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Research  on  Deterministic  Methods  of  Seismic  Source  Identification 

Summary  of  Semi-Annual  Technical  Report  for 
October  1,  1982  -  Sept.  30,  1984 

The  objectives  of  the  research  conducted  during  the  2  year  contract  were  to:  (1)  Develop 
and  test  methods  of  discrimination  in  the  regional  and  teleseismic  distance  range  using  physical 
source  parameter  discriminants,  (2)  Pursue  theoretical  and  observational  studies  of  seismic 
sources;  (3)  To  develop  methods  of  theoretical  seismogram  synthesis  in  the  near,  regional  and 
teleseismic  distance  ranges  for  structure  and  source  definition;  (4)  To  develop  and  apply 
advanced  signal  processing/analysis  methods  for  discrimination  and  explosion  yield  estimation 
studies  and;  (5)  to  pursue  near  field  studies  of  explosions  and  earthquakes  for  detailed  source 
definition. 

In  this  report  we  describe  specific  research  results  pertaining  to:  (1)  The  theoretical  basis 
for  automatic  seismic  signal  detection  and  analysis,  and  (2)  Analytical  methods  for  the  represen¬ 
tation  of  seismic  radiation  fields  in  uniformly  layered  elastic/anelastic  media.  This  modal 
method  provides  predictions  of  both  body  and  surface  waves  in  the  frequency  range  from  0  to 
about  15  HZ  at  near  and  regional  distances  from  seismic  sources.  This  latter  exposition  is 
intended  to  be  comprehensive  and  integrates  new  and  old  results  and  methods.  The  modal 
representation  method  for  seismic  radiation  fields  is  being  employed  to  describe  radiation  fields 
for  the  prediction  of  earthquake  and  explosion  radiation  fields  and  has  been  used  to  evaluate  a 
variety  of  detection  and  discrimination  methods. 

Previous  annual  and  semi-annual  reports  have  described  applications  of  the  QHD  signal 
analysis  methods  and  the  seismic  synthesis  methods  to  source  discrimination  problems.  In  par¬ 
ticular,  these  earlier  reports  have  described  research  results  from  the  observational  studies 
related  to  discrimination,  wherein  computer  programs  employing  the  theoretical  results 
described  here  are  used  to  analyze  the  data  and  provide  the  basis  of  interpretation  of  the  data 
for  regional  discrimination  and  yield  estimation. 
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Part  1  -  Methods  of  Seismic  Signal  Detection 
and  Analysis  Based  on  Quasi-Harmonic 
Decomposition  (QHD) 


I.  INTRODUCTION 


Any  time  series  can.  of  course,  be  decomposed  into  Fourier  or  harmonic 
components,  and  in  this  case  the  outcome  is  that  we  determine  the  amplitude 
and  relative  phase  of  each  harmonic  component  occuring  within  the  time  win¬ 
dow  chosen  for  such  a  decomposition.  We  do  not,  however,  know  anything  about 
when  the  energy,  associated  with  a  particular  frequency,  may  be  arriving  within 
this  time  window.  That  is.  we  have  no  time  resolution  within  the  window.  If  we 
attempt  to  sub-divide  the  time  window  and  perform  a  Fourier  analysis  of  these 
subsets  of  the  time  series  in  order  to  obtain  a  higher  degree  of  time  resolution 
then,  of  course,  truncation  effects  begin  to  severely  contaminate  the  amplitude 
spectral  estimates. 

It  is.  however,  well  known  that  information  about  the  time  of  energy  arrival 
(i.e.,  the  group  time)  can  be  retained  if  one  is  prepared  to  accept  a  certain 
amount  of  uncertainty  in  the  frequency  value  to  be  associated  with  the  spectral 
amplitude  (or  energy  magnitude)  and  phase  estimate  (e.g.,  Helstrom.  1960).  In 
particular,  simple  narrow  band  filtering  of  the  time  series,  x(t).  at  a  set  of  filter 
center  frequencies  /*,  can  be  used  to  generate  an  associated  set  of  quasi¬ 
harmonic  time  series,  y(t\f  t ),  having  properties  such  that  the  spectral  ampli¬ 
tudes  and  phases  of  pulses  making  up  the  original  time  series  can  be  estimated, 
along  with  the  tine  at  which  each  pulse  spectral  component  has  arrived  within 
the  time  window.  That  is,  the  spectral  character  and  arrival  time  of  each  fre¬ 
quency  component  of  the  individual  "signal"  or  "noise"  pulses,  making  up  the 
original  time  series,  may  be  estimated. 

The  price  to  be  paid  for  this  simultaneous  estimate  of  both  frequency  con¬ 
tent  and  (energy)  arrival  time  is  an  uncertainty  in  the  precise  time  and  fre¬ 
quency  to  be  assigned  to  the  spectral  components  of  each  pulse  or  transient 
making  up  the  original  time  series.  Thus  if  hf  is  the  bandwidth  of  one 


particular  narrow  band  filter  at  its  half  power  point,  then  the  uncertainty  in 
the  estimate  of  the  energy  arrival  time  is  related  to  A/  by  an  “uncertainty  prin¬ 
ciple"  in  the  form.*: 

—  —  1 

At)  Af  — 

where  A t  is  the  uncertainty  (variance)  in  the  estimate  of  the  energy  or  group 
arrival  time  and  At)  —  2nA /.  (e.g.  If  a  Gaussian  narrow  band  filter  is  employed, 
then  the  uncertainty  product  is  a  minimum,  so  that  At)  A t  =1/2.  For  all 
other  niters  At)  At  >  1/  2.) 

In  the  context  of  this  “quasi-harmonic  decomposition"  of  the  time  series  by 
narrow  band  filtering,  we  may  think  of  Fourier  decomposition  as  the  limiting 
case,  in  which  the  narrow  band  filters  that  are  convolved  with  the  time  series 
have  zero  bandwidth.  Thus  Fourier  decomposition  corresponds  to  the  case 
At)  =  0  in  (1),  and  in  this  limiting  case.  At  becomes  unbounded.  This  simply 
means  that  we  obtain  precise  frequency  information  and  no  information  about 
energy  arrival  time,  within  the  time  window  chosen  for  analysis. 

The  possibility  of  time  resolution  as  well  as  frequency  resolution  in  analyz¬ 
ing  a  time  series,  however,  allows  “instantaneous"  measures  of  the  spectral 
character  of  the  series,  that  is  measurements  of  time  varying  phase  and  spec¬ 
tral  amplitude  as  well  as  measurements  of  derived  quantities  such  as  the 
“instantaneous  frequency",  corresponding  to  the  time  derivative  of  the  time 
varying  phase.  These  quantities  have  been  discussed  extensively  in  the  litera¬ 
ture  (e.g.  Helstrom.  1980)  and  are  meaningful  when  defined  for  quasi-harmonic 
time  series.  Thus,  our  objective  here  is  to  show  how  decomposition  of  a  broad¬ 
band  time  series  into  a  set  of  quasi-harmonic  time  series,  by  Gaussian  narrow 

*  8m  Denny  and  Chin  (1976),  for  example,  for  a  more  preciM  definition  of  the  quantities 
At)  and  At  appearing  in  (1). 


band  filtering,  can  be  optimally  accomplished  so  as  to  allow  appropriately  accu¬ 
rate  measurements  of  time  varying  phase,  amplitude  and  frequency  for  each 
quasi-harmonic  component  to  be  made.  We  will  then  show  how  this  derived  set 
of  time  dependent  variables  may  be  used  to  detect  "desired  signals",  with 
prescribed  physical  properties,  in  the  presence  of  "undesired  signals"  or 
"noise".  Specifically  we  will  show  how  a  variety  of  "signal  sensitive"  filters  can  be 
defined  to  isolate,  for  example,  signals  with  particualr  dispersive  properties 
(dispersion  filters),  wave  number  and  frequency  characteristics  (  o—k  filters) 
and  polarization  or  particle  motion  characteristics  (polarization  filters). 
Further,  we  can  use  these  instantaneous  variables  to  define  another  class  of 
"filters"  designed  to  separate  interfering  (time  overlapping)  signals,  with  these 
latter  interference  detection  filters  being  based  on  the  variance  of  the  instan¬ 
taneous  frequency  variable  and  the  splitting  (jump  discontinuities)  in  the 
dispersion  characteristics  of  a  multiple  pulse  signal. 

The  joint  use  of  all  such  filters  proves  to  be  a  very  powerful  "matched 
filter"  for  the  isolation  or  decomposition  of  complex,  multiple  pulse,  time  series. 
In  addition  to  signal  detection  however,  the  decomposition  of  the  time  series 
into  quasi-harmonic  components  and  the  generation  of  time  varying  spectral 
variables  allows  the  whole  time  series  to  be  analyzed  in  considerable  detail.  In 
particular  signal  spectra,  with  noise  corrections,  can  be  obtained  along  with 
frequency  dependent  signal  dispersion,  polarization  and  wave  number  vector 
direction  and  magnitude.  These  are  essentially  all  the  signal  properties  that  we 
need  to  (or  can)  know  in  order  to  interpret  the  time  series  in  some  physical 
context.  In  this  regard,  we  will  show  how  deterministic,  as  well  as  probablistic 
noise  corrections,  may  be  applied  to  isolated  signed  data,  and  how  this 
corrected  signal  data  may  be  used  to  make  physical  inferences  regarding  the 
origins  of  the  signals.  To  do  this  we  need  to  consider  a  specific  kind  or  type  of 


time  series,  in  order  to  not  only  relate  it  to  its  physical  origins,  but  also  in 
order  to  define  the  kinds  of  "signals"  to  be  detected. 

Thus,  while  the  general  signal  analysis  methods  discussed  may  be  applied 
to  any  type  of  time  series  regardless  of  its  origins,  we  will  focus  on  seismic  wave 
generated  time  series  since,  aside  from  being  familiar  to  us.  they  are  particu¬ 
larly  interesting  in  view  of  their  uses  (e.g.  seismic  methods  for  resource 
exploration,  event  identification  for  underground  nuclear  event  detection, 
planetary  structure  investigations,  earthquake  prediction  studies)  and  since 
they  are  particularly  rich  in  structure  (a  vector  wave  field,  at  least  two  intrinsic 
wave  propagation  velocities,  and  a  variety  of  guided  wave  possibilities  in  plane¬ 
tary  structures  leading  to  complexities  in  signal  dispersion,  polarization  and 
wave  number  vector  direction  and  magnitude).  In  the  application  of  QHD 
methods  to  the  seismic  data  discussed  in  this  study,  we  will  consider  examples 
of  automated  signal  detection  and  analysis  that  relate  directly  to  event  identifi¬ 
cation  and  earth  structure  determinations. 

Numerous  previous  studies  have  addressed  the  analysis  of  seismic  data 
using  related  multiple  narrow  band  filtering  techniques.  In  particular,  in  an 
early  study,  Alexander  (1963)  used  multiple  band  pass  filtering  to  determine 
the  group  velocity  dispersion  of  seismic  surface  waves,  from  which  elastic  velo¬ 
city  structure  in  the  earth  was  inferred.  Arcbambeau  ef  a l  (1965),  Archambeau 
and  Flinn  (1965)  and  Archambeau,  Flinn  and  Lambert  (1966)  used  multiple  nar¬ 
row  band  filtering  to  determine  the  dispersion  and  spectral  properties  of  both 
body  and  surface  waves  for  earth  structure  and  source  studies.  Dziewonski, 
Block  and  Landisman  (1969)  applied  multiple  filtering  extensively  to  give 
seismic  surface  wave  dispersion  measurements,  and  studied  the  resolution  of 
the  method  using  synthetic  seismograms.  Dziewonski,  Mills  and  Block  (1972) 
Investigated  the  effect  of  preprocessing  strongly  dispersed  surface  wave  signals 
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with  a  matched  filter  in  the  time  domain,  followed  by  narrow-band  filtering,  and 
showed  that  this  procedure,  which  they  termed  a  residual  dispersion  measure¬ 
ment,  was  superior  to  the  direct  measurement  of  the  total  dispersion.  Several 
later  studies,  for  example  by  Inston,  Marshall  and  Blamay  (1971),  Cara  (1973) 
and  Denny  and  Chin  (1976),  have  also  considered  schemes  for  optimization  of 
the  multiple  filtering  method  in  order  to  achieve  greater  accuracy  and  resolu¬ 
tion  in  seismic  surface  wave  dispersion  measurements. 

Most  of  this  work  has  focused  on  signal  dispersion  measurements,  most 
commonly  for  low  mode  order  seismic  surface  waves.  In  the  present  study  we 
will  generate  signal  dispersion  results  in  a  manner  similar  to  that  used  in  these 
earlier  studies,  but  with  an  optimization  approach  that  invokes  matched  filter¬ 
ing,  when  required,  along  with  Alter  bandwidth  variations  with  frequency.  We 
will  then  go  a  step  further  and  specify  dispersion  filters  whose  function  is  to 
isolate  signals  with  particular  dispersion  (or  residual  dispersion)  characteris¬ 
tics.  In  doing  this  we  will  also  adjoin  to  the  dispersion  filter,  other  filters,  whose 
function  is  to  simultaneously  search  for  other  required  signal  properties,  in 
particular  polarization  and  wave  vector  directional  properties. 


In  regard  to  seismic  polarization  filters,  or  particle  motion  detectors,  early 
work  focused  on  time  domain  analysis  without  the  explicit  use  of  special  multi¬ 
ple  narrow  band  filtering.  Examples  of  such  filters  are  those  discussed  by 
Shimshoni  and  Smith  (1964)  and  Mims  and  Sax  (1965).  Flinn  (1965)  used  time 
domain  particle  motion  measurements  to  perform  both  polarization  and  direc¬ 
tional,  or  wave  number  vector,  filtering.  Arcbambeau  and  Flinn  (1965)  defined 
polarization  filters  in  the  spectral  domain  for  seismic  body  and  surface  wave 
detection  and  analysis.  In  later  work,  Archambeau,  Flinn  and  Lambert  (1969) 
applied  these  polarization  filters  as  a  means  of  isolating  multiple  body  wave 
arrivals  for  study  of  the  earth's  upper  mantle  structure.  Body  wave  isolation 
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and  enhancement,  with  suppression  of  surface  wave  noise  by  polarization  tech¬ 
niques,  has  been  subsequently  used  by  a  number  of  investigators,  for  example 
by  Montalbetti  and  Kanasewich  (1970),  who  used  an  extension  of  the  time 
domain  approach  described  by  Fiinn  (1965);  and  by  Lewis  and  Meyer  (1968).  who 
employed  the  frequency  domain  approach.  On  the  other  hand,  Simons  (1968) 
applied  frequency  domain  polarization  filtering,  using  a  moving  window  FFT,  to 
detect  and  enhance  seismic  surface  waves.  (These  earlier  polarization  filtering 
methods,  and  their  applications,  are  summarized  in  more  detail  by  Kanasewich 
(1975).)  Further,  both  Fiinn  (1965)  and  Simons  (1968)  considered  directionality 
filters  using  the  amplitude  of  the  ground  motion  recorded  on  the  three  spatial 
components  defining  the  vector  wave  field  at  a  single  three  component  seismic 
station.  This  filtering  amounts  to  "vector  wave  number  filtering"  at  a  single 
point,  since  the  wave  number  vector  describes  the  normal  direction  to  the  wave 
front  and,  for  a  prescribed  wave  type,  the  recorded  components  of  ground 
motion  can  be  used  to  determine  the  orientation  of  this  vector.  Thus,  only  cer¬ 
tain  desired  wave  number  vector  orientations  are  accepted  by  such  a  direc¬ 
tionality  filter. 

In  the  present,  study  we  will  combine  polarization  filtering  with  wave 
number  vector  filtering.  Polarization  filters  can  be  defined  in  terms  of  the 
instantaneous  phase  variable  generated  by  quasi-harmonic  decomposition 
(QHD)  of  the  individual  vector  components  of  a  seismic  wave.  In  this  case  the 
phase  difference  between  two  spatial  components  of  the  wave  field  serves  to 
define  the  polarization.  Thus  we  can  obtain  a  frequency  dependent,  time- 
varying  measure  of  the  polarization  of  the  seismic  wave  signals  and  noise  mak¬ 
ing  up  the  time  series.  Similarily,  we  can  define  vector  wave  number  directions 
in  terms  of  the  time  varying  amplitude  components,  taken  at  the  group  arrival 
times;  where  we  obtain  the  spectral  information  at  the  group  arrival  times  and 
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define  directionality  that  is  both  frequency  and  time  dependent.  Compared  to 
previous  time  domain  methods,  this  approach  has  the  advantage  of  defining  a 
polarization  filter,  and  wave  number  vector  filter,  for  each  filter  center  fre¬ 
quency,  so  that  frequency  dependent  polarization  and  directionality  filtering 
can  be  performed.  In  addition,  greater  time  resolution  can  be  achieved  by  the 
QHD  method  than  is  possible  by  the  previously  used  spectral  methods,  and  trun¬ 
cation  effects  arising  from  moving  spectral  windows,  are  entirely  avoided. 
Further,  we  can  combine  dispersion  measurements  with  polarization  and  direc¬ 
tionality  measurements  in  one  single  operation,  with  the  QHD  procedure  gen¬ 
erating  all  the  required  information.  Then  optimization  of  the  filter  design  for 
any  one  of  the  measurements,  for  example  for  accuracy  in  signal  dispersion 
measurements,  ensures  that  the  filter  design  is  also  optimal  for  all  the  other 
measurements,  such  as  for  time  varying  amplitude,  phase,  polarization,  etc. 
Most  important,  however,  is  the  fact  that  we  may  associate  a  set  of  signals  (or 
noise)  related  variables,  such  as  spectral  amplitude,  phase,  polarization  and 
instantaneous  frequency,  with  each  measured  group  arrival  time,  which  itself 
corresponds  to  the  time  at  which  seismic  wave  energy,  at  (or  near)  the  filter 
center  frequency,  arrives  at  the  sensor.  This  generation  of  time  connected 
functions  of  frequency,  which  may  be  associated  one  to  one  with  signal  (or 
noise)  pulses  in  the  time  series,  provides  us  with  the  opportunity  of  performing 
a  series  of  joint  filtering  operations.  These  filtering  operations  are  both  time 
varying  and  frequency  dependent  in  general,  and  may  be  linear  or  non-linear, 
but  are,  in  any  case,  used  to  extract  desired  signals  along  with  estimates  of 
their  spectral  properties.  In  the  next  section  we  will  consider  the  appropriate 
seismic  signal  definitions  so  as  to  establish  a  physical  context  in  which  these 
filtering  operations  can  be  defined,  implemented,  and  applied. 
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n.  SEISMIC  SIGNAL  DEFINITIONS 

In  the  previous  discussion  we  indicated  that  "joint  filtering"  operations  are 
to  be  considered.  What  is  meant  is  that  a  number  of  criteria,  based  on  known 
physical  characteristics  of  the  various  types  of  seismic  signals,  will  be  used 
together,  to  detect,  or  search  for,  signals  of  a  desired  type.  This  is  a  matched 
filter  concept  with  signals  detected  on  the  basis  of  correlations  with  specified 
signal  properties.  However,  our  application  of  the  matched  filtering  approach 
will  be  quite  different  than  the  usual  linear  filtering  method,  which  involves 
generation  of  what  amounts  to  a  cross-correlation  of  the  specified  signal  wave 
form  with  a  time  series,  and  selection  of  correlation  peaks  above  a  threshold  as 
signal  detections.  Instead,  we  will  define  only  general  "robust”  characteristics 
of  the  signal  types  of  interest,  rather  than  attempt  to  specify  details  of  wave 
forms  or  spectra,  and  then  use  quasi-harmonic  decomposition  methods  to  (con¬ 
tinuously)  monitor  the  time  series  for  the  signal  characteristics  of  interest. 
Detections  of  "desired  signals"  would  then  be  obtained  at  those  time  intervals 
during  which  the  time  series  exhibited  all  (or  most  of)  the  characteristics  of  the 
specified  signal  type.  Obviously  how  closely  and  in  what  sense  the  time  series  is 
to  match  the  discrete  set  of  specific  signal  properties  must  also  be  precisely 
defined  --  and  this  will  be  considered  in  detail  in  a  later  section.  At  this  point 
our  task  is  to  specify  signal  definitions,  and  then  to  devise  the  filter  analysis 
proceedures  which  will  provide  data  for  detections.  Once  a  desired  signal  form 
is  detected,  the  niter  output  data  can  also  provide  details  of  the  signal  spectra, 
timing  and  so  on. 

The  signal  characteristis  that  will  be  used  as  criteria  for  detection  and  iso¬ 
lation  of  seismic  signals  are  those  invariant  properties  that  are  usually  used  to 
define  the  various  seismic  body  and  surface  wave  types.  Specifically  we  have; 
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(A)  Seismic  body  wave  pulses  characterized  by: 

(1)  Minimal  dispersion;  such  that  each  frequency  component  of  the 
signal  pulse  has  nearly  the  same  group  and  phase  velocity. 

(2)  Linear  polarization;  with  compressional  body  waves  (P  waves)  hav¬ 
ing  particle  motion  in  the  direction  of  the  wave  number  vector 
and  normal  to  the  wave  front;  and  with  shear  waves  (SV  and  SH) 
having  particle  motion  in  directions  perpendicular  to  the  wave 
number  vector,  in  the  plane  tangent  to  the  wave  front. 

(3)  Pulse-like  amplitude  spectrum;  where  the  amplitude  spectrum 
has  broad  band  width  and  is  a  slowly  varying  function  of  fre¬ 
quency. 

(4)  Minimal  variation  in  wave  number  orientation;  with  nearly  con¬ 
stant  direction  (azimuth  and  emergence  angles)  of  the  wave 
number  vector  as  a  function  of  frequency. 

(5)  Near  vertical  wave  vector  directions;  so  that  the  apparent  emer¬ 
gence  angle  is  relatively  low  and  the  apparent  phase  velocity  is 
large. 

(B)  Seismic  surface  waves,  characterized  by; 

(1)  Strong  dispersion;  such  that  the  group  velocity  generally 
increases  with  decreasing  frequency,  due  to  the  generally 
increasing  elastic  velocities  with  depth  in  the  earth. 

(2)  Elliptical  polarization  ( Rayleigh  type  surface  waves)  or 
transverse  linear  polarization  (Love  type  surface  waves);  with 
fundamental  mode  Rayleigh  waves  usually  being  retrograde  ellipt- 
ically  polarized.4,  so  that  the  particle  motion  is  in  elliptical  orbits, 
with  the  onset  of  motion  along  the  wave  number  axis  in  the  oppo¬ 
site  direction  from  the  wave  number  vector.  Higher  mode  Rayleigh 
waves  may  be  either  retrograde  or  prograde  elliptically  polarized, 
the  latter  with  elliptical  particle  motion  initiating  in  the  same 
direction  as  the  wave  number  vector. 

For  all  Rayleigh  modes,  the  eccentricity  of  the  particle  orbits  —  or 
the  ellipticity  of  the  polarization  —  are  functions  of  frequency, 
having  a  variation  which  is  similar  in  form  to  that  for  the  wave 
group  velocity. 

(3)  Peaked  amplitude  spectrum:  so  that  the  time  domain  signal  asso¬ 
ciated  with  a  particular  surface  wave  mode  is  usually  well 
dispersed  and  extended  in  time,  with  at  least  several  cycles  of 

4 Fundament*]  mod*  Rayleigh  waves  can,  in  certain  structures,  within  particular  frequency 
bands,  also  be  prograde  elliptically  polarised. 
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motion  present.  Thus  surface  waves  typically  tend  to  be  narrower 
in  band  width  than  body  waves  and  appear  to  "ring"  in  the  time 
domain. 

(4)  Moderate  variation  in  wave  number  azimuth  orientation  -with  fre¬ 
quency;  with  higher  modes  showing  more  variability  than  funda¬ 
mental  mode  surface  waves,  these  variations  being  principally 
associated  with  lateral  variations  in  the  near  surface  velocity 
structure  within  the  earth. 

(5)  Horizontal  wave  number  vector  directions;  so  that  phase  velocities 
are  low  compared  to  those  for  body  waves. 

These  characteristics  are  the  major  fixed  properties  of  seismic  "phases" 
and  so  are  the  most  important  criteria  for  phase  detection  and  identification. 
The  signal  properties  are,  out  of  necessity,  rather  broadly  defined  but,  as  we 
will  show,  they  are  sufficient  to  provide  the  basis  for  automated 
detection/identification  of  signals  in  a  high  noise  background.  Specifically,  the 
approach  is  to  specify,  as  precisely  as  possible,  what  type  of  signal  is  to  be 
detected  in  terms  of  wave  type  (body  wave;  Rayleigh  wave,  etc.)  and  the  dis¬ 
tance  and  azimuth  range  for  the  originating  event,  so  that  signals  with  the 
appropriate  dispersion,  polarization,  wave  number  vector  orientation  and  spec¬ 
tral  properties  can  be  isolated.  The  first  filtering  process  to  be  applied  is  simply 
one  designed  to  continuously  generate,  as  functions  of  time,  the  required  fre¬ 
quency  dependent  dispersion,  polarization,  apparent  wave  vector  and  spectral 
information;  with  the  detection  operation  then  being  a  selection  process, 
wherein  one  picks  out  the  signals  having  desired  properties.  It  is  our  purpose  to 
make  this  selection  process  an  objective,  computer  controlled,  operation. 
Naturally  the  details  of  this  process  can  be  quite  involved  and  will  depend,  at 
least  in  part,  on  the  nature  of  the  noise  background.  It  is  therefore  necessary 
to  also  consider  the  character  of  the  expected  noise,  or  undesired  signal,  in 
order  to  effectively  design  signal  isolation  procedures. 
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m.  NOISE  DEFINITIONS:  SEISMIC 


In  most  applications  it  is  desirable  to  be  able  to  detect  and  isolate  particu¬ 
lar  types  of  signals  from  seismic  sources  within  certain  prescribed  distance 
ranges  and  azimuth  sectors.  For  example,  it  is  often  most  important  to  be  able 
to  detect  compessional  (P)  waves  from  sources  at  teleseismic  distances,  occur- 
ing  within  some  azimuth  sector,  and  to  reject  signals  from  nearby  microearth¬ 
quakes  as  well  as  teleseismic  signals  from  events  outside  the  azimuth  sector  of 
interest.  Further,  it  is  important  to  be  able  to  reject  all  background  seismic 
noise  due  to  more  or  less  randomly  distributed  sources,  such  as  are  associated 
with  atmospheric  pressure  fluctuations,  ocean  surf,  and  cultural  effects. 
Clearly  undesired  signals  from  isolated  seismic  sources  outside  the  distance 
and/or  azimuth  range  of  interest  may  only  differ  from  the  signals  of  interest  in 
wave  number  vector  orientation,  or,  if  near  the  receiver,  in  high  frequency  con¬ 
tent.  On  the  other  hand,  characeristics  of  the  background  seismic  noise,  while 
obviously  having  wave  characteristics  that  are  similar  to  those  previously 
enumerated  for  "signals",  do  have  details  of  spectral  and  modal  structure  that 
are  quite  distinctive.  These  properties  provide  the  basis  for  filtering  operations 
which  can  distinguish  between  most  signals  of  interest  and  background  noise. 

The  noise  characteristics  can  be  described  as  follows: 

(1)  General  spectral  properties  and  modal  composition:  The  average  back¬ 
ground  noise  spectrum,  obtained  from  long  time  noise  samples,  is 
strongly  peaked  in  the  period  range  from  6  to  8  seconds,  with  a  secon¬ 
dary  peak  of  much  lower  level  near  0.3  seconds.  Other  secondary  low 
level  maxima  may  occur  at  very  long  periods,  but  roughly  speaking, 
the  noise  at  lower  frequencies  decreases  from  the  peak  at  6  -  8 
seconds,  rapidly  at  first,  then  much  more  slowly.  At  high  frequencies, 
f>10  Hz,  the  noise  is  mostly  wind  generated  and  originates  locally,  with 
very  high  levels  possible.  Essentially  all  the  low  frequency  noise, 
including  that  associated  with  the  sharp  spectral  peaks  at  8  -  8 
seconds  corresponds  to  fundamental  mode  Rayleigh  type  surface 
waves,  with  minor  contributions  from  higher  modes  and  from  Love 
type  surface  waves.  The  noise  peak  near  0.3  seconds  corresponds 
mainly  to  higher  mode  Rayleigh  waves,  with  most  of  the  noise  peaks 
being  composed  of  a  superposition  of  many  higher  modes  and 
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therefore  having  body  wave-like  characteristics.  Thus  at  these  higher 
frequencies,  and  particularly  near  the  0.3  sec.  spectral  peak,  the  noise 
is  predominantly  higher  mode  Rayleigh  type  surface  waves,  with  some 
body  waves  contributing  as  well. 

(2)  Spatial  and  temporal  variability  in  noise  spectral  levels,  invariance  of 
mode  composition:  The  overall  level  of  the  noise  spectrum  is  dependent 
on  proximity  to  noise  sources,  such  as  ocean  boundaries,  and  is  also 
quite  variable  with  time.  However,  the  modal  compostion  of  the  noise  is 
effectively  constant,  even  though  the  overall  mode  excitation  level 
may  be  quite  variable.  (This  implies  that  all  the  surface  wave  modes, 
fundamental  and  higher,  are  excited  by  a  typical  noise  source  and 
that  most  of  the  fundamental  mode  excitation  is  at  6  to  8  sec.  and 
most  higher  mode  excitiation  is  near  0.3  sec.) 

(3)  Pulse-like  time  domain  composition  of  high  frequency  background 
noise:  Analysis  of  the  detailed  time  domain  structure  of  the  back¬ 
ground  noise,  for  example  by  narrow  band  filtering,  indicates  that, 
especially  in  the  intermediate  to  high  frequency  range,  the  noise  can 
be  described  as  a  superposition  of  discrete  noise  "pulses"  consisting  of 
a  few  cycles  of  motion.  In  view  of  the  normally  random  space-time  dis¬ 
tribution  of  sources,  a  noise  time  series  is  therefore  composed  of  a 
random  time  sequence  of  "bursts"  or  pulses,  with  each  pulse  having  a 
mode  characteristic  of  the  type  previously  described.  Each  pulse  is 
therefore  actually  a  coherent  signal  propagating  as  a  surface  wave  or 
occasionally  as  a  body  wave. 

The  essential  differences,  between  the  noise  and  the  signals  we  wish  to 
detect,  are  that  the  body  wave  signals  will  be  undispersed,  linearly  polarized 
pulses  with  wave  number  vectors  that  are  oriented  within  a  narrow  spatial 
"cone"  and  with  associated  high  values  of  apparent  phase  velocity.  On  the  other 
hand,  surface  wave  noise  has,  for  the  most  part,  moderate  to  strong  dispersion, 
is  elliptically  polarized,  has  narrow  band  width  compared  to  the  signals,  and  has 
relatively  low  apparent  phase  velocity.  Thus  filtering,  designed  to  exploit  the 
differences  in  dispersion,  polarization,  frequency  content  and  wave  number 
characteristics,  will  be  effective  in  isolating  body  wave  signals.  Further,  that 
part  of  the  noise  field  having  body  wave  characteristics  will  have  random  wave 
number  vector  orientations  and  can  be  eliminated,  in  pari,  by  wave  number 
filtering. 

On  the  other  hand,  surface  wave  signals  of  interest  are  typically  of  broader 
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band  width  than  the  surface  wave  noise,  which  is  strongly  peaked  in  the  6-8 
sec.  period  range,  and  so  ordinary  frequency  filtering  can  be  used  to  reduce 
the  surface  wave  noise  without  appreciable  loss  of  signal  power.  The  noise  also 
has  random  wave  number  vector  orientation  characteristics,  and  both  vector 
wave  number  filters  and  matched  filters  can,  therefore,  be  used  to  select  sur¬ 
face  wave  signals  from  particular  spatial  directions. 

IV.  ANALYTICAL  REPRESENTATIONS  FOR  SEISMIC  SIGNALS  AND  NOISE 

The  signals  and  the  time  series  we  are  dealing  with  are  associated  with  vec¬ 
tor  wave  fields  and  so  are  themselves  vector  quantities.  Ve  will,  in  general,  deal 
with  one  spatial  component  of  the  field  at  a  time  and  will  suppress  component 
indices  when  the  treatment  applies  equally  to  all  signal  components.  Ve  will 
only  be  explicit  as  to  vector  components  when  it  is  necesary  to  combine  them, 
as  is  the  case,  for  example,  when  defining  single  receiver  polarization  or  wave 
number  filters. 

Modal  Representation  of  Pulse-like  Signals  and  Noise 

Any  component  of  an  observed  seismic  displacement  time  series,  x(t),  can 
be  expressed  as  sum  of  "signals",  in  the  form 


*(0  =  £  5,(0  =  £  Re/Jf  &(«)*<“*<£/  (l) 

II  II 

where  Re  denotes  the  real  part  of  the  Fourier  integral.  Here,  since  x(t)  and 
Sn(t)  are  real,  then  with  the  Fourier  transform  of  5^(0  denoted  by  Sti(o).  we 
have  3*(-«)  =  3£(<y).  where  is  the  complex  conjugate  of  3*-  Similarity, 

MS  Ms  MS 

*  (— «),  =  x  (&>)•  where  x(o)  is  the  Fourier  transform  of  x(t). 

The  sum  in  (l)  is  over  "signals"  S*(t)  which  may,  in  fact,  be  noise  pulses 


depending  on  how  we  choose  to  define  what  we  call  a  signal.  At  this  point  we  can 
simply  call  an  energy  pulse  in  the  time  series  a  signal  of  some  sort.  This  is  con¬ 
sistent  with  our  previous  descriptions  of  seismic  signals  and  noise,  where  both 
were  described  as  having  specific  wave  type  properties.  We  can,  as  we  will 
demonstrate  below,  represent  each  "pulse"  5*,(f)  as  a  superposition  of  modes, 
whether  it  is  noise  or  signal  and  whether  it  corresponds  to  a  body  or  surface 
wave.  If  SJj  happens  to  correspond  to  a  body  wave,  then  its  representation  is  a 
superposition  of  many  higher  modes.  If  it  is  a  surface  wave,  then,  of  course,  it  is 
represented  by  just  one  mode.  In  any  case  the  "signals"  all  have  the  general 
form  of  propagating  vector  fields  and,  in  the  mode  representation  form,  are 
given  in  the  frequency  domain  by: 

&(«)  =  E  4P  (r*«)  e*P  -*[*«»■  +sC(«)]J  (2) 

where  the  sum  over  m  is  a  sum  over  modes  and  r  is  the  earth’s  radial  coordi¬ 
nate.  Here  AJ1  is  the  complex  wave  number  for  the  mode  of  the  n &  signal. 
The  phase  p™  can  be  called  the  initial  phase  for  the  mode  term  of  the  n— 
signal,  and  is  associated  with  the  phase  delay  at  the  source.  The  function  is 
source  and  medium  dependent,  but  because  of  its  strong  dependence  on  the 
source  of  the  signal,  will  be  called  a  mode  excitation  function.  The  mode  excita¬ 
tion  function  for  the  n&  signal,  A£\  may  of  course  be  zero  (or  negligible)  for  all 
except  one  mode  index  value,  so  that  the  (infinite)  sum  over  m  can  be  degen¬ 
erate.  For  body  waves  this  is  not  the  case  however  and,  for  a  given  body  wave 
signal,  many  of  the  Aj?  functions  will  be  of  comparable  magnitude. 

The  detailed  nature  of  body  wave  modal  structure  is  illustrated  by  the 
mode  curves  of  Figure  (l)  and  (2).  from  Harvey  (1980).  Seismograms  generated 
from  the  excitation  of  these  modes  by  a  theoretical  earthquake  source  are 
shown  in  Figure  (3). 
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Figure  1.  A  section  of  the  modal  phase  velocity  versus  frequency 
curves  for  the  velocity-density  structure  in  the 
Southern  California  region,  just  east  of  the  San 
Andreas  fault  zone  and  south  of  Los  Angeles  (structure 
after  Kanamori  and  Hadley,  1978) .  The  prominent 
flattening  of  the  higher  mode  curves  near  4.6  km/sec, 
again  in  the  range  around  6.2  to  6.5  km/sec  and  near 
8.0  km/sec  produce,  respectively,  S,  P«  and  body 
wave  arrivals  when  the  modes  are  summed,  as  shown  in 
Pigure  3. 


Numerous  body  wave  pulses  are  generated  by  the  superposition  of  these 
modes,  but  each  seismic  body  phase  is  associated  with  a  large  number  of  higher 
modes  all  having  a  flat  segment  in  their  phase  velocity  dispersion  curves  at 
nearly  the  same  phase  velocity  value.  This  is  illustrated  in  Figure  1  by  the  flat 
dispersion  in  the  modes  near  the  phase  velocity  of  8  km/sec.  (This  corresponds 
to  Pn,  the  crust-mantle  boundary  refracted  wave  which  is  the  first  arrival  on 
the  theoretical  seismogram  shown  in  Figure  3.)  When  the  phase  velocity  is  con¬ 
stant,  as  it  is  on  the  flat  parts  of  the  dispersion  curves  in  Figure  1,  then  the 
phase  and  group  (energy)  velocities  are  equal,  and  both  constant.  Figure  2 
shows  the  more  complex  group  velocity  curves,  where  their  constant  group 
velocity  near  8  km/sec.  is,  nevertheless,  also  evident.  Thus,  for  the  superposed 
modes,  the  flat  sections  of  the  phase  velocity  curves  will  contribute,  in  sum.  to  a 
pulse  like  seismic  phase,  with  all  frequency  components  arriving  at  nearly  the 
same  time  (ie.,  with  nearly  constant  group  velocity).  Because  mode  excitation  at 
a  given  frequency  is  inversely  proportional  to  the  derivative  of  the  group  velo¬ 
city  with  respect  to  frequency,  then  when  the  group  velocity  is  nearly  constant, 
as  it  is  along  the  flat  parts  of  the  dispersion  curves,  its  derivative  is  small  and 
the  excitation  of  th.e  mode  will  be  high.  Since  there  are  many  modes  with  the 
same  locally  constant  group  and  phase  velocities,  then  not  only  will  they  all  sum 
to  give  an  undispersed  pulse,  but  each  mode  will  have  relatively  larger  excita¬ 
tion  in  this  (narrow)  group  velocity  range  than  it  does  at  either  higher  or  lower 
values  of  group  velocity,  where  the  dispersion  curves  are  steep.  Thus,  one  can 
expect  a  well  defined  pulse,  with  larger  amplitude  than  that  of  the  more 
dispersed  energy  arriving  just  before  and  after  it  in  time. 

This  description  constitutes  a  modal  representation  of  a  classical  body 
wave.  In  the  other  extreme,  the  theoretical  seismograms  in  Figure  (3)  show  late 
arriving  energy  corresponding  to  the  single  fundamental  mode  Rayleigh  wave. 
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Figure  2.  Low  frequency  sections  of  the  phase  and  group  velocity  curves  as 
functions  of  frequency.  Only  group  velocity  dispersion  curves  for  the 
fundamental  and  first  few  higher  modes  are  shown,  due  to  the  complexity  of 
these  functions  in  the  velocity-frequency  plane.  Maxima  and  minima  along  the 
modal  group  velocity  curves  give  rise  to  "Airy  phases",  where  the  mode 
excitation  is  relatively  much  larger  than  at  other  points  on  the  curves, 
where  the  slope  is  larger. 
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Early  tine  section  of  the  synthetic  seismograms 
generated  by  mode  superposition  for  the  Southern 
California  structure  shown  in  Figure  1.  The  lines 
are  theoretical  travel  times  for  the  "body  wave" 
phases  Pn  and  Pgr  based  on  simple  ray  theory.  The 
complex  P„  wave  train  results  from  many  higher  mode 
contributions  in  the  phase  velocity  range  around 
6.2  to  6.5  km/sec ,  where  the  mode  curves  show  low 
slopes  and  flattening  as  functions  of  frequency  as 
illustrated  in  Figure  1.  The  sharp,  time  compressed, 
pulse  results  from  the  superposition  of  the  large 
number  of  modes  with  very  pronounced  flattening  in 
a  narrow  range  near  8.0  km/sec. 


bo  that  for  the  representation  of  this  signal,  only  a  single  mode  is  involved. 

Much  of  the  seismic  energy  propagation  must,  however,  be  described  as 
intermediate  between  a  classical  body  wave,  composed  of  a  very  large  number 
of  higher  modes  with  near  constant  phase  and  group  velocities,  and  ordinary 
fundamental  mode  surface  waves.  An  example  of  this  is  the  Pg  wave  energy, 
arriving  in  the  group  velocity  range  from  around  6.2  to  5.0  km/sec  in  the  syn¬ 
thetic  seismograms  shown  in  Figure  3.  The  origins  of  this  large  amplitude,  time 
distributed  arrival  are  higher  mode  contributions,  associated  with  the  relatively 
slight  mode  curve  flattening  in  the  phase  velocity  range  from  about  7  to  9 
km/sec.  In  this  case  there  is  no  extremely  well  defined  coincidence  of  flat  sec¬ 
tions  in  the  mode  curves  at  a  constant  phase  and  group  velocity  (giving  a  well 
defined  body  wave  that  can  be  approximated  by  a  ray)  but  instead  one  finds  a 
less  pronounced  flattening  which  occurs  over  a  range  of  phase  and  group  velo¬ 
cities.  In  this  case  the  engergy  transmission  cannot  be  viewed  very  simply  in 
terms  of  rays,  as  are  classical  body  phases,  but  must  be  viewed  as  either  a  large 
number  of  higher  modes,  contributing  over  a  group  and  phase  velocity  range  or 
as  a  large  number  of  rays  with  steep  angles  of  incidence  (high  phase  velocities) 
representing  multiple  reflections  and  refractions  distributed  in  arrival  time 
over  a  relatively  long  time  interval. 

Mode  Groups  as  Isolated  Signals 

A  modal  description  is  a  formal  analytical  representation  that  will  be  used 
here  to  provide  a  mathematically  complete  representation  of  the  entire  seismic 
wave  field.  We  will,  when  appropriate  however,  make  use  of  the  special  features 
of  the  modal  superposition  when  dealing  with  body  wave  signals.  Now,  for  such  a 
representation,  we  can  combine  (1)  and  (2)  to  obtain: 


(3) 


x{t)  -  £  Re  X.  iC  (r.w)exp 


9k, m 
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where  is  the  phase  factor  for  ER.  Thus,  from  previous  definitions,  every 
"signal"  Sn(t )  is  represented  by  a  sum  of  modal  contributions  Ef£  such  that 


&  (r.o)  =  2  (*•“)  (5_b) 

m 

Here,  as  elsewhere,  m  is  restricted  to  a  set  of  integers  (Zn)  appropriate  to  the 
ntk  signal.  Precisely  how  we  chose  to  define  the  signals  Sn(t)  will  determine 
which  modes,  and  which  phase  or  group  velocity  segments  of  the  modes,  contri¬ 
bute  to  a  given  signal. 

The  wave  number  AJI*  is  conveniently  represented  in  terms  of  a  phase  velo¬ 
city  C£*  for  the  m—  mode  of  the  n—  signal,  so  that 


AT(«)  *  6 (6) 

In  this  representation  of  an  observed  seismic  time  series  we  will  use  a  rela¬ 
tion  like  (3)  to  describe  the  entire  time  series  and  to  represent  a  superposition 
of  "signals"  from  many  sources,  including  sources  which  produce  "signals"  that 
are  not  of  interest,  as  well  as  signals  from  some  particular  source  which  are  of 
interest.  Thus  (3).  and  equivalent  mode  sums,  represent  a  multiple  source 
superposition,  with  the  sum  over  n  representing  a  sum  over  different  source 
signals  as  well  as  a  sum  over  different  signals  from  the  same  source.  In  this  case 
only  a  few  signals  from  the  series  are  important,  while  all  other  terms  in  the 

l 

i 


sum  over  n  are  viewed  as  noise  relative  to  the  signals  to  be  isolated.  Our  objec¬ 
tive  is  to  devise  systematic  methods  that  may  be  used  to  pick  out  of  this  series 
those  terms  that  satisfy  prescribed  conditions,  such  as  those  for  body  waves  of 
a  particular  type  from  a  source  in  a  particular  distance-azimuth  range. 

V.  GAUSSIAN  NARROW  BAND  FILTERING  OF  SEISMIC  TIME  SERIES 

The  essential  tool  to  be  used  to  isolate  desired  signals  from  the  seismic 
time  series  will  be  the  quasi-harmonic  decomposition  of  the  time  series.  This 
decomposition  provides  time  varying  spectral  information  such  that  particular 
characteristics  of  the  sought  for  signals  can  be  measured  and  used  to  automat¬ 
ically  select,  or  filter  out,  the  appropriate  signals. 

Narrow  Band  Filter  Combs 

In  order  to  decompose  the  time  series  we  will  use  a  "comb",  or  set.  of  ultra 

narrow  band  Gaussian  filters,  with  center  frequencies  fk,  k=l,2 . N.  The 

transfer  function  of  a  Gaussian  filter  is  defined  to  be: 

a  <««)  *  Vf  (7) 

where  ok  =  2nfk  is  the  angular  center  frequency  of  the  filter,  and  a  is  a 
parameter  determining  the  bandwidth  of  the  filter.  The  impulse  response  of  the 
filter  is  given  by 

8k(t)  -  rl  [ft  it)  =  ^  /_7  [>/ “ 

where  F~l  denotes  the  inverse  Fourier  Transform  operator.  We  therefore  have, 
after  evaluating  the  integral: 

*(0  =  <B) 

as  the  Impulse  response.  Hence  this  is  a  simple  sinusoid  modulated  by  an 


exponential  function  of  time,  in  both  time  directions  from  t  =  0.  Thus  this  filter 
is  phaseless  and  when  applied  to  filter  a  signal,  it  will  introduce  no  phase  shift 
and  will  "ring"  in  both  the  positive  and  negative  time  directions.  Clearly  this 
modulation  can  be  controlled  by  the  choice  of  a.  (Although  not  explicitly  indi¬ 
cated,  a  can  be  chosed  to  be  a  function  of  v*  for  the  set  of  filters  ( Gu  ).  ) 


A  useful  means  of  specifying  the  filter  bandwidth  is  by  use  of  the  filter 
quality  factor,  or  Q,  which  is  defined  to  be  the  ratio  /*  /A/* .  with  A/*  equal  to 
half  the  bandwidth  at  the  half  power  point  of  the  filter.  To  express  the  parame¬ 
ters  a  in  terms  of  the  filter  Q,  we  observe  first  that  the  power  at  u  =  ok  is  given 
by: 


ato  =  a* , 


so  that  the  half  power  points,  at  u  —  fi>*+Ao *,  are  such  that 


|c4(w*  ±A»fc;a)|*  =  |^j-j 

Evaluating  explicitly  the  expression  for  G»,  at  o  -  gives  the  condition 

that: 


In  terms  of  the  half  bandwidth  Av*.  we  have  for  a: 


Using  the  definition  for  the  filter  Q,  then: 


Thus  the  set  of  Gaussian  filters  have  the  form 


(9-a) 


(9-b) 


(10-a) 
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Analytical  Forms  of  Gaussian  Filter  Responses 
for  Isolated  Signals 

If  we  apply  any  one  of  the  Gaussian  filters  to  a  signal  S^(t ),  then  the  out¬ 
put  of  the  filter  can  be  denoted  as  (t;ok),  and  is  given  by  the  inverse  Fourier 
transform  of  the  product  of  the  signal  transform  with  the  filter  transfer  func¬ 
tion,  that  is: 

5.  (*:&>*)  =  &  (u)  fnt(»-6>*;Qt)e<ulif» 

Here  we  have  written  out  the  functional  dependence  of  the  filter  on  the  differ¬ 
ence  o  —ok,  and  on  the  parameter  Q*.  explicitly.  In  order  to  evaluate  this 
integral  we  introduce  a  change  of  variable  such  that 


o’  =  o-ok 


so  we  get 


<*>*»  4-  1 

5*  («;«*)•=  Re  -gjj- 3*  («'+«*)  («’;&)  e^diy’j  (11) 

and  then  consider  the  expansion  of  the  signal  spectrum  §t(u)  about  the  fre¬ 
quency  o  -  Of  that  is  about  the  point  o'  -  0. 

We  observe  that  we  have  from  (5): 

3,  («)  =  £  A?  (r,o)e~*  *^rM) 

m 

with  the  sum  running  over  the  modes.  Thus,  for  any  signal  S*  we  can  expand 
and  i 'J*  individually  about  o  -  ok,  and  obtain: 

rfTfr.u’+B,)  *  £  it  [sS'jC<r’")l,.<"'>' 
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*n(r .«•+«*)  sp  Ji  *r(*’.o)J^(w')1 

We  will  see  that  if  the  bandwidth  of  the  filters  is  narrow  enough,  then  only  small 
values  of  o'  contribute  to  Sn(t, ot).  and  the  expansion  for  S r,  resulting  from  the 
expansions  in  (12).  will  converge  rapidly. 

In  order  to  simplify  the  notation  somewhat,  we  note  that  the  mode  sum 
over  m  commutes  with  the  integral  over  tS  in  (11.).  Thus  the  filter,  in  effect, 
acts  on  each  mode  independently  and  the  total  output  is  just  the  sum  of  the 
individually  filtered  modes.  For  this  reason  we  will  temporarily  suppress  the 
mode  index  m  and  the  mode  sum.  since  the  manipulations  and  integral  evalua¬ 
tion  are  independent  of  which  mode  is  being  considered.  In  this  case  we  can.  for 


example,  write 


*il)  (r.uk)  =  jj-  (r-u)]ot 


as  a  notational  convenience  in  the  second  expansion  in  (12.) 
We  observe  from  (5.)  that  (with  the  index  m  suppressed): 


Therefore, 


+»  (*■.«*)  *  K*  r  +  fn(w) 


Since  the  group  velocity  Un(u)  is  given  by 


then 


*.*.**>.>  V  .•»  .> 
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*i0>  (r,uk)  =  K*  (wk)  +  fn(«k)  = 


&>k 


♦  ?«(«■>*) 


where  C*(uk)  is  the  phase  velocity  evaluated  at  u~uk  ■ 

Introducing  group  and  phase  delay  times,  tg  and  tp,  defined  by 


*in)M  =  Un(ok)  *  y"0(a>t) 
W-W  ,  r  .  ?*(«*) 

=  cST  + 


(15) 


uk 


for  the  signal  (with  the  mode  index  again  suppressed)  and  then  using  the 
expansions  of  (12)  in  the  integral  relation  (11)  for  the  filter  output  S*(t;ok), 
gives: 


5»(*:«*)  =  Re 


i_  rj  -V'-’i 
Ar  «  2tt 
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■  £  rr  ATOM* )  7  («')'  -  «’("  ) 
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exp 
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where 


W  -  £ 


P«8 


cf  w* 


To  evaluate  the  integrals  in  this  expansion  for  £*  we  can  express  the  factor 
exp  \-iP(u)\  in  a  power  series  in  u  and  use  the  integral  relation: 


/(*»)"  « w  ■to*  =  V7  kf/H :  Re(o)  >  0 


to  give  (16)  in  series  form. 


The  exponential  expansion  required  has  the  form 


E(o)  m  exp[- */»(«)]  =  1+1)  •««’ 

m«2 


where  the  first  few  coefficients  ^  are: 

d*t}n 
do2 
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Thus,  the  filtered  output  in  (16)  can  be  expressed  as: 


=  Re 
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or.  combining  the  two  summations,  as: 


=  R« 
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£tf0)(r.o*)  =  i4j0,(r.wk)  ;  A»{I)(r.o*)  =  i4n(,)(r.«*) 
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:  P*  2 


Now  using  the  integral  relation  previously  quoted,  we  get,  after  grouping 
terms  in  order  by  powers  of  Lok : 

*«W*)  =  [^-jRe  Je<tt*{‘  [^°)(r.at) 

♦  2 [^^(r.o*)  ♦  <(/  -  t0(n))  B&l)(r.ok)  ] 

+  4  (  3#W(r.c*)  +  3i<f  -  t}n>)  ^(r.ur)  -  (t  -  tjn')z  Bp\ r.ok  )  Au*4 

+  0(A«*#)  ]  exp  -  ~~  (t  -  tj"))*  ]  }  (18) 

Here,  in  terms  of  the  fundamental  signal  parameters,  the  coefficients  B&p)  are: 

^0)(r.o*)  =  ^(r.Ofc) 


BPH r.ok)  =  ~  Al*Hr,ok)  -  Ai9)(r,uk ) 

B&*>(r.»k)  =  'APW)  -  |-j^0)(r.«t)  (^5-^ 


♦  3A,(1)(r 


£^4)(r.«k)  a  ~  [^(r.o*)  -3^,(r.«*)  ] 


-26- 


RSI 


V  {f  Si 

-  E  «~l  -* 


■->  K.  I 


+  8 ApK 


where  the  superscript  index  on  the  signal  amplitude  An  denotes  the  order  of 
derivative  with  respect  to  angular  frequency,  that  is: 


AilKr.uk) 


l  do1  ]«* 


This  expression  for  the  Filtered  signal  is  formally  exact,  under  the  tacit  assump¬ 
tion  that  the  expansions  of  the  signal  amplitude  and  phase  functions  are  con¬ 
vergent.  Such  convergence  is.  however,  assured  for  a  physically  realizable  sig- 


There  are  a  number  of  important  features  in  equation  (16.)  that  deserve 
comment.  First  of  all.  the  result  is  similar  to  the  impulse  response  of  the  filter 

itself,  in  that  it  consists  of  a  sinusoid,  namely  modulated  by  a  decay¬ 

ing  exponential  in  time  -  that  is  modulated  by 


In  addition,  however,  the  result  contains  a  second  complex  valued  modulation 
function,  involving  the  signal  spectrum  and  frequency  derivatives  of  the  spec¬ 
trum.  plus  frequency  derivatives  of  the  group  delay,  t£n\  all  evaluated  at  the 
filter  center  frequency  .  This  function  therefore  carries  information  con¬ 
cerning  the  signal  spectrum  and  its  dispersive  characteristics,  whereas  the 
pure  sinusoid  carrier  involves  only  the  phase  delay  of  the  signal  and  the  Gaus¬ 
sian  exponential  modulation  depends  only  on  the  group  delay  of  the  signal. 
Further,  this  function  has  the  form  of  a  power  series  in  so  that  for  A &>* 
small  it  can  be  approximated  by  the  first  few  terms.  It  is  therefore  useful  to 


'r 


write  the  filtered  signal  result  in  the  form: 


where 


S*(t.uk)  -  (  jrr)Re  ^(t -tp>;Aok)  •  **  (< ** 
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Thus,  the  modulation  function  14*  is  a  relatively  single  power  series  in  Aw*  with 
coefficients  which  are  polynomials  in  the  time  difference  (f  —  It  is  also 

useful  to  express  in  polar  form,  in  which  case  we  can  write  the  filter  output 

as: 

§»(#:w*)  -  _^n)),cos[«*(<  -*f]  (19a) 
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(20a) 


=  tan-1 

A +  0(  A«jf) 

Here,  in  expressing  these  last  results,  we  have  assumed  A}P  A  0. 


(20b) 


limiting  Forms  of  the  Quasi-Harmonic  Filtered  Signal.  5*  - 

An  important  property  of  the  modulation  function  is  that 

lim  Ac*)  =  Ai0)  (rc*)  (21) 

Jim  5,(<*>.)  =  (^-)  R «  (22) 

so  that,  as  the  filter  bandwith  approaches  zero,  we  obtain  a  pure  sinusoid  with 
amplitude  and  phase  equal  to  the  Fourier  amplitude  and  phase.  This  of  course 
is  not  unexpected,  since  this  limit  must  correspond  to  the  operation  of  ordinary 
Fourier  analysis.  Nevertheless  this  shows  that  we  may  approach  the  limit  of 
Fourier  spectral  analysis  if  Ac*  is  made  small,  and  we  will  of  course  wish  to 
obtain  good  spectral  estimates  of  the  signal  by  narrow  band  filtering,  as  well  as 
accurate  group  arrival  time  estimates,  by  choosing  Ac*  small.  In  this  regard 
however,  the  limiting  case  given  in  (22.)  shows  that  we  can  obtain  exact  spectral 
estimates  by  narrow  band  filtering  when  A&*-»0,  corresponding  to  Fourier 
decomposition,  but  no  estimate  of  the  group  arrival  time. 

We  observe  however,  from  (20.),  that  when  t  =f/n)  for  the  n&  signal,  then 


most  of  the  terms  in  the  expression  vanish,  and: 
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A!8)  dt  )  _ 

Un  (0;Aq*)  =  1  +  ^jor  -  i  Ad*8  +  O(Loi)  (23) 

with  the  final  terms  being,  at  most,  of  order  Lot  and  involving  products  of  high 
order  derivatives  of  An  and  with  respect  to  frequency.  Therefore  when  A u* 
is  such  that 

A  uf[<«p(n)]  ,  % 

7Fl'^rK<<1  (24) 

then  we  can  expand  y  in  powers  of  this  small  factor,  which  gives 
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for  the  7  factors  appearing  in  the  function  Afn.  These  factors  can  be  written  in 
polar  form  as 

- i _ 


7~»  w  e  2  ' 


7~2  « 


(26-a) 
(26- b) 


where 


We  see  that  as  long  as  the  condition  (24)  holds,  then  the  factors  y~l  and  7"® 
appearing  in  Afn  can  be  approximated  as  simple  phase  shifts.  In  this  situation, 
the  value  of  at  t  =t^  is,  from  (23): 

Consequently  the  narrow  band  Gaussian  filtered  signal  S*,.  at  has  the 

form: 
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(27) 

for  Ap)?0  at  Comparing  this  to  (22),  the  limiting  case  of  the  Fourier 

decomposition  for  which  Ao*  is  zero,  we  now  observe  that  we  get  a  very  similar 
result  at  t  =tp)  when  Aw*  is  finite,  so  long  as: 
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In  this  case  we  note  that  the  phase  factor  $g  in  (27)  can  be  neglected 
since. 
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by  the  first  inequality  in  (28).  Further  when  the  second  inequality  in  (28)  also 
holds,  the  complex  amplitude  factors  in  (27)  can  be  neglected  relative  to  unity. 
Thus  we  can  clearly  approximate  the  limiting  case  in  which  the  amplitude  of  the 
filtered  output  is  equal  to  the  Fourier  amplitude  at  the  center  frequency  of  the 
niter  when  tok  is  finite,  if  we  measure  the  amplitue  of  the  filtered  time  series  at 
t=tp)  ,  the  group  arrival  time.  Similarity  the  phase  of  the  filtered  time  series  at 
t=tp)  is,  to  a  good  approximation,  equal  to  ok  (tp)-tp))  which  is  the  Fourier 
phase  of  the  n—  signal  at  the  frequency  u*. 


Filter  Design  and  Signal  Conditioning  Requirements 


In  view  of  the  previous  limiting  case  results,  we  observe  that  the  objective 
of  narrow-band  filtering,  which  is  to  decompose  a  complex  time  series  into  a 
series  of  nearly  harmonic  components  that  simultaneously  retain  accurate  sig¬ 
nal  energy  arrival  time  and  Fourier  amplitude  and  phase  information  for  each 
signal,  can  be  achieved  provided  the  filter  design  and  signal  characteristics  are 
such  that  the  Q  of  each  filter  simultaneously  satisfies  the  conditions 


Qk » |42,(o* )/  240)(t>* )  |,/2 


(29) 


where  Q  and  ok  are  the  quality  factor  and  center  frequency  for  the  filter. 
For  the  filter  design  criteria  to  be  met,  the  expected  signals  within  the  time 
series  should  be  minimally  dispersed  and  have  smooth  amplitude  spectra.  Body 
wave  signal  pulses  intrinsically  have  such  characteristics,  while  other  seismic 
"phases"  require  preprocessing  of  the  time  series  to  render  them  of  the 
required  form.  Specifically,  preprocessing  with  a  matched  filter,  which  is  the 
approximate  inverse  of  the  expected  signal,  within  the  bandwidth  to  be  covered 
by  the  set  of  narrow  band  filters,  is  required. 

An  acceptable  design  of  the  narrow  band  filter  set  is  achieved  using: 

Qk  -rbok 

where  r*  is  a  constant  such  that: 

(^T)m‘  ;  \^2)M/ZA^\ok)\  } 

In  this  case  the  bandwidth,  &ok.  is  constant  for  all  the  filters  of  the  set.  since 
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from  the  definition  of  the  filter  Q  we  have 

buk  - 1/  Tb 

Hence,  r*  is  a  time  like  parameter  which  is  a  measure  of  the  ring  time  of  the 
narrow  band  filters,  (ie.  The  amplitude  of  the  impulse  response  of  each  filter  is 
decreased  by  1/e,  from  its  peak  at  when  t=tg± Tb.  )  In  addition,  the 

resolution  time,  Af*.  for  any  time  measurement  from  a  filter  output  is  con¬ 
strained  by  the  "uncertainty  principal",  having  the  form;  bok  Af*  >_t,  for  all  k. 
Here  Ao*  is  again,  the  filter  bandwidth  and  c  is  some  constant  independent  of 
the  filter  index  k.  Consequently  the  uncertainty  in  measuring  the  time  at  which 
a  filter  output  has  a  maximum  is  proportional  to  rk.  in  particular: 

A tk  2.  trb 

Since  we  cannot  precisely  define  e  ,  other  than  to  say  that  it  is  greater  than 
zero  and  at  most  of  order  1,  then  we  cannot  say  much  more  than  that  Af* 
increases  or  decreases  directly  with  t*.  With  the  choice  of  the  filter  set  Q  fac¬ 
tors  as  indicated  above  however,  then  every  Filter  of  the  set  has  the  same 
theoretical  resolution,  as  well  as  bandwidth.  That  is,  A tk  and  Aw*  are  both  con¬ 
stant  for  all  k  and  directly  or  inversely  proportional,  respectively,  to  r* . 


niter  Output  For  a  Time  Series  of  Superposed 
Pulse-like  Signals  and  Noise 

When  the  narrow  band  filter  design  meets  the  conditions  specified  by  (29), 
then  the  filter  output,  for  a  input  series  of  signals  £y,(f  ),  has  the  form* 


KmU  -<*(w,;Ao*)  *  AjSi(ok)  +  U&Ko*  )(t  -te™)+%A™(ok  )]  +  oJao*]  (31) 

Here  the  mode  sum  over  m  has  been  restored  to  the  expressions  for  complete¬ 
ness.  The  mode  index  m  takes  on  values  from  sets  (J*)  ,  which  are  determined 
by  the  signal  index  n.  (  ie.  Different  sets  of  modes  are  summed  for  different  sig¬ 
nals.)  The  function  Y(t;uk)  is  introduced  to  denote  the  k—  narrow  band  filter 
output  of  a  time  series  consisting  of  a  sum  of  signals  of  pulse-like  form. 

The  results  expressed  in  (30)  and  (31)  therefore  apply  to  signals  of  pulse 
form,  that  is  minimally  dispersed  and  spectrally  smooth.  In  general  the  time 
series  is  composed  of  signals  of  this  type,  plus  noise,  the  latter  being  all  "sig¬ 
nals"  not  of  this  form.  Thus,  in  general,  we  have  a  time  series  x(t)  which  is  a  sum 
of  signals  5^(0  and  noise  Ni(t)  .  where  the  noise  is  indexed  to  indicate  that  it 
too  is  composed  of  a  series  of  discrete  wave  packets  (typically  overlapping  in 
time)  which  are  usually  dispersed  and/or  have  "non-smooth"  spectral  ampli¬ 
tude  character.  Nevertheless  the  noise  can  be  represented  by  modes,  however 
these  modes  are  not  excited  by  the  single  source  responsible  for  the  signals, 
but  rather  by  a  more  or  less  random  distribution  of  sources  scattered  in  space 
and  time.  Thus  we  write  the  basic  time  series  as: 


*(0  =  £5»(0  +  5>t(f) 

n  I 

and  the  output  of  the  narrow  band  Gaussian  filter  as 

X(t;ot)  =  £^(f;«*)  +  £*,(*;«*)  (32) 

n  I 

where  the  First  sum  on  the  right  side  is  Y(t  ;ok).  as  given  in  (30)  -  (31).  and  the 
second  sum,  representing  the  filtered  noise,  can  be  denoted  Z(t;o *).  From  the 
previous  general  analysis  it  is  quite  evident  that  Z(t;vk)  can  be  expressed  in 
the  same  general  form  as  is  Y(t  \ok ).  that  is  as 


Z(t;uk)  m  :t>k)  =  2^7 Re 


2 
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Here,  however,  the  modulation  function  now  has  a  complex  form,  given  for¬ 
mally  by  equation  (20),  which  cannot  in  general  be  reduced  to  the  simple  form 
of  (31).  since  the  various  frequency  derivatives  of  the  group  time  and  spectral 
amplitudes  for  the  noise  generally  do  not  conform  to  the  smoothness  properties 
required.  Further  the  group  and  phase  times,  tg  and  tp  .  are  generally  different 
for  each  mode  contributing  to  a  particular  noise  packet  or  wave  group.  Hence 
these  times  are  explicitly  indexed  by  both  1  the  "wave  packet"  index,  and  m  the 
index  for  the  modes  associated  with  that  wave  packet. 

We  note  however,  that  both  the  group  time  t and  the  phase  time  tjn^ ,  for 
the  n—  signal  pulse  appearing  in  Y(t  ;&>*),  are  nearly  constant  and  equal  to  each 
other.  In  particular,  each  mode  that  contibutes  to  a  given  signal  pulse  will  have 
values  of  ta  and  tp  that  are  constant  over  the  frequency  range  for  which  the 
mode  has  significant  excitation.  Since  each  contributing  mode  will  have  essen¬ 
tially  the  same  values  of  ta  and  tp  for  a  particular  signal  pulse,  then,  for  seismic 
body  wave  signals.  ta  and  tp  are  mode  independent,  and  so  are  not  indexed  by  m 
in  (30)  -  (31).  The  same  observations  hold  for  wave  types  that  are  not  pulse-like 
initially,  but  that  have  been  rendered  pulse-like  and  undispersed  by  matched 
filtering.  In  this  case-the  dispersed  modes  of  interest  have  been  mapped,  by  the 
filtering,  into  nearly  undispersed  pulse-like  arrivals. 

Consequently,  for  such  signals  we  can  formally  carry  out  the  sum  over 
modes  indicated  in  (30),  taking  account  of  the  fact  that  the  exponential  terms 
do  not  depend  on  m  for  seismic  body  waves  and/or  matched  filtered  signals, 
and  rewrite  the  results  in  the  simpler  form: 

r{t  wk)*j^Re  £  "n (t  -t}n)*uk  )« '  (< V°>  (34) 


with 


(35) 


Here  the  amplitude  factors  A^X4*)*  representing  the  p&  order  frequency 
derivatives  of  the  amplitude  spectrum  of  the  signal  pulse,  correspond  to 
sums  over  only  the  modes  contributng  to  the  n—  signal.  That  is: 


-£*&(<*)  (36) 

m 

where  m  is  restricted  to  the  mode  index  set  ( L *)  that  corresponds  to  the  n— 
signal.  (We  do  not,  a-priori,  know  what  mode  set  does  in  fact  belong  to  a  partic¬ 
ular  signal,  but  it  is  sufficient  for  the  present  to  know  that  such  finite  sets 
exist.  At  a  given  frequency  it  is  typical  that  only  one  or  two  modes  contribute  to 
a  particular  body  wave  phase  -  as  illlustrated  in  Figure  1.  Thus,  the  sum  in  (36) 
is  generally  over  only  one  or  two  mode  indices,  m.  at  a  fixed  frequency  v>t  ■  ) 


The  principal  use  of  the  narrow  band  Filtering  output  is  to  provide  meas¬ 
urements  of  signal  group  arrival  time  (  i,  values)  and  the  amplitudes  and  phase 
of  the  signals  of  these  times.  Thus  the  essential  feature  of  the  output  from  a 
set  of  narrow  band  filters  is  the  behaviour  of  the  filter  output  at,  or  near,  the 
energy  or  group  arrival  times  ta.  At  a  group  time  of  a  particular  signal,  say  the 
j—  signal,  a  properly  designed  filter  with  center  frequency  has  an  output: 


ArWW) «  i- 


COS  tik(tjjto—t^ fn>) 


♦  Re  s{jh>( V” - )« .....]  (37) 

1 

where  i Mim  in  the  final  sum  is  the  modulation  function  for  the  noise,  and  has  the 
form  given  by  (20). 

The  first  term  in  this  expression  is  just  equal  to  the  Fourier  component  of 


the  signal  at  tt*  ,  while  the  second  term  is  a  sum  over  all  other  (n?j  )  pulse 
like  signals,  and  the  final  term  is  the  sum  over  all  the  noise  arrivals.  The  latter 
two  terms  can  be  viewed  as  interference  terms  relative  to  the  signed,  whose 
properties  (  eg.  spectrum)  are  to  be  estimated.  Clearly,  because  of  the  Gaus¬ 
sian  exponential  factor  contained  in  both  these  terms,  the  sums  will  be  small 
relative  to  the  first  term  if 

(38) 

In  this  case 

and  we  get  an  uncontaminated  result  for  the  j—  signal.  If  on  the  other  hand, 
the  time  separation  between  this  signal  and  other  signals  and  noise  is  not  very 
large,  then  some  of  the  terms  in  the  two  series  may  contribute  significantly  to 
the  value  of  X  at  t  -  t ^  .  However,  measurements  of  X  at  a  signal  group 
arrival  time  can  kc  corrected  for  such  "contamination"  to  give  an  estimate  of 
the  Fourier  component  of  the  signals,  and  equation  (37)  can  be  used  as  the 
basis  for  this  correction.  The  details  will  be  considered  in  a  later  section. 

VI.  MEASUREMENTS  OF  FILTER  OUTPUT  CHARACTERISTICS  USING  ASSOCIATED 
FILTER  FUNCTIONS:  ENVELOPE.  INSTANTANEOUS  PHASE  AND  FREQUENCY  FUNCTIONS 

In  order  to  obtain  signal  spectrum  and  dispersion  estimates  from  the  nar¬ 
row  band  filtering  operations,  it  is  convenient  to  define  a  set  of  "associated" 
filter  functions,  derived  from  the  analytical  forms  obtained  for  the  filter  output, 
which  permit  simple  computer  controlled  determinations  of  signal  and  noise 
group  arrival  times  and  amplitude  spectra.  Clearly  the  most  important  task, 
considering  the  emphasis  placed  on  measuring  the  filter  output  at  group  times 
tg,  is  to  be  able  to  determine  the  t0  values  for  all  signal  and  noise  arrivals 
within  the  time  series.  Ve  will  therefore  first  consider  the  filter  envelope 


function,  from  which  we  can  obtain  estimates  of  the  t9  values. 

The  Quadrature  Signal 

As  a  preliminary,  it  is  necessary  to  define  the  "quadrature  signal",  which 
can  be  formed  from  the  original  time  series  spectrum,  as 

£(«)  =  — i  sgn  ( u )  X(i>)  (39) 

where  sgn(o)  denotes  the  algebraic  sign  of  o.  and  X(o)  is  the  Fourier  spectrum 
of  the  original  time  series,  while  £(t»)  denotes  the  spectrum  of  the  quadrature 
time  series.  In  the  time  domain,  x(t)  is  the  Hilbert  transform  of  z(t).  The  pro¬ 
perty  of  x  which  is  useful  here  is  that  its  Fourier  components  are  shifted  in 
phase  by  — n/2.  as  is  easily  seen  from  (39).  This  is  the  origin  of  its  designation 
as  the  "quadrature  signal". 

How,  for  any  time  series,  *(f  )•' 

*<0  =  ^  f^X(u)*iu,dU  =  f’  [*(«)««* ^(-^e^Jdu 

However,  since  x(t)  is  real,  then  X(u)  must  have  the  property  that 
X{— «)  =  Jf* (u),  where  X*  («)  is  the  complex  conjugate  of  X(u).  Thus,  with 

X(o)  =  |*(»)|e«w. 

then  it  follows  from  this  property  that: 

x(t)  =  “To  |*(w)|cos[w*  +  *  (w)]d«  =  —ReJ^  |[A'(a>)  eiutdo  (40) 

Similarily  the  quadrature  time  function  x(t),  with  spectrum  given  in  (39). 
is  represented  by: 

But  again,  since  X{—u)  =  X*(u),  we  have 
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or 

x  =  ~^f0  |*(w)|cos[of  +  t-n/2)do  -  ^-R efo  |jA'(t>)|e<{*(")""/2)je<"‘dcj  (41) 

(Comparison  of  this  Fourier  expansion  for  x(t)  with  that  for  x(t )  now  explicitly 
shows  that  x(t)  has  Fourier  components  that  are  shifted  in  phase,  from  those 
for  x(t ).  by  a  factor  of  — n/Z.) 

The  Complex  or  Analytic  Signal 

The  quadrature  signal  can  be  used,  together  with  the  original  signal,  to 
define  what  is  usually  termed  the  complex  or  analytic  signal  .  A  complex  signal 
representation  provides  the  basis  for  the  definition  of  instantaneous  phase  and 
frequency  functions,  as  well  as  an  envelope  function.  Specifically,  we  define  the 
complex  signal.  xe(t ).  as: 

*«(f)-*(0  +i*(f)  (42) 

and  from  this  definition  it  follows  immediately  .  from  the  representations  (40) 
and  (41)  for  x(t)  and x(t),  that 

*c(0=  /,"*(«)•“* do  —  \*i{ut  **]d o  (43) 

Clearly.  xe(<  )  is  a  complex  valued  function  of  time,  rather  than  being  real  as  are 
sr(f  )  and  x(t),  and.  in  particular,  can  be  written  in  the  forms 

*e(f)*Rej*e(0]  +  i  Im[*c(f  )]=  |xe(<)|e<#e(<)  (44) 

where 

•,  «tan“l 


(45) 


Thus.  *,(0.  can  be  written  in  polar  form,  since  it  is  complex,  and  can  easily  be 
computed  using  (43).  which  is  simply  the  transform  of  x(t)  over  only  the  posi¬ 
tive  frequencies.  Further,  it  follows  from  these  definitions  that 


*  (O  =  Re 


ee(oj  =  |*e(0|cos[*£(0] 
r(f)  =  Im|xe(<)j=  |*e(0|  sin  [mo] 


(46) 


Thus,  an  arbitrary  time  series,  that  may  not  have  sinusoid  form,  can  neverthe¬ 
less  be  expressed  in  this  Nsinusoid  form"  through  the  use  of  its  quadrature 
function. 

In  analogy  with  the  ordinary  frequency  of  a  sinusoid,  the  phase  factor, 
4>c(f).  appearing  in  (44)  -  (46),  can  be  expressed  in  terms  of  some  fixed  fre¬ 
quency,  say  oq.  and  a  frequency  modulation  function,  +( t ).  through  the  defini¬ 
tion: 


(47) 


Thus,  the  time  derivitaive  of  fe  corresponds  to  a  frequency,  in  particular  the 
instantaneous  frequency  0(f),  and  we  have  from  (47) 

0(0=*  a,0  =c*,  +  *«)  (48) 

In  terms  of  the  "carrier"  frequency  &>0  and  frequency  modulation  function  't(t). 

From  the  definiton  of  4e(0  in  (45)  we  also  have 


mo* 


(49) 


Envelope,  Instantaneous  Phase  and  Frequency  Functions 
For  Quasi-Harmonic  Time  Series 


For  a  quasi-harmonic  time  series  the  amplitude  and  phase  functions  in  (45) 


are  usefully  defined  as  the  envelope  and  instantaneous  phase  of  the  quasi¬ 
harmonic  time  series.  That  is,  with  x*(t;ofc)  denoting  the  quasi-harmonic  time 
series,  the  envelope  and  instantaneous  phase  are  defined 


Ek(t;uk)s  j*f  (<;«*)]  = 

[*t  (*;«•*) 


(50) 


)=?£(*:«*)  =  tan 


-i 


**(<  ;«*) 


where  |af  |  and  are  the  amplitude  and  phase  of  the  complex  signal  formed, 

using  (43).  from  the  Gaussian  filtered  time  series  z(t).  Similarity  the  instan¬ 
taneous  frequency  associated  with  each  quasi-harmonic  time  series  xk(t :ok)  is 


0* (*:«*)  = 


d<pt 


dt 


dxk  ^  dxk 


Xk~dT~Xk  dt 


/  E? 


For  a  time  series  consisting  of  pulse-like  signals,  then  we  have: 


zk(t;ok)  =  £$,(*;«*) 


where,  from  (20) 


(*:«*)  *  ^ Re 


-tuf 


Slu*  ^  Ji*  «  -tM;£x>k )  e  '**+  {‘ 
with  the  phase  time  variable,  defined  in  (15),  having  the  form 


•^mskr* 


*«(«*) 


ok 


and  with  the  modulation  function,  given  by  (20),  having  the  form 


(51) 


(52) 


(53) 


-  41  - 


(54) 


Here  tbe  factor  y  is  such  that: 


where 


•#  =  tan-1 


A 

~W 


as  shown  in  (24)- (28).  Now  we  can  express  in  polar  form  as 


*.  =  |n,  |  v*i 


(55) 


where 


^  |  =  >u£0) 


1  4 


A**> 

Aof 

2i4n(0) 

20 

4  0(Awif) 


fn  *  tan-1 


W 


+0(5ufl 


and  rewrite  3*  as: 


(56) 


$»(*;«*)  = 


IK  I  (i -*;»>)* 


2n 


.  .  f.i  ^ a 

COSCJ*  4 - — S- 

*  w*  2uk 


(57) 


Similarly,  the  quadrature  signal  has  the  form 


t>t  Zuk 


(58) 
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All  these  expressions  are  accurate  up  to  terms  of  order  We  observe,  how¬ 
ever,  that  the  phase  of  Un,  as  expressed  by  above,  is  only  meaningful  when  t 
is  near  that  is  when  |f  — <1.  It  is  also  true  that  the  higher  order  terms 
in  IM.I  involve  (f  - 1 }n)),  and  these  terms  become  large  as  t  - 1 becomes  large. 
Thus  while  the  Gaussian  exponential  factors  in  S*,  and  produce  a  rapid  con¬ 
vergence  of  these  functions  to  zero  when  t  - 1 ^  is  large,  it  is  nevertheless  clear 
that  the  analytical  expressions  for  these  functions,  involving  truncated  power 
series  in  (t  are  only  accurate  for  t  near  tjn\ 

We  now  observe  that  the  envelope  function,  £*(f  ;&>*),  for  the  whole  filtered 
time  series,  can  be  related  analytically  to  envelope  functions  for  the  individual 
“signals".  In  particular,  (see  also  equation  (  )): 

=  2  (59-,) 

Vi 

with 


H(n>(t;uk)  m  |55|  =  S* 


(59- b) 


and 


xtH*'**)  *  ?«(<:<*)  =  tan'1 


(59-c) 


where  |<Sgj  is  the  modulus  of  the  “complex  signal"  associated  with  the  filtered 
signal  pulse  $»,  *nd  **  the  phase. 

Now  from  (57)  and  (58)  we  heve. 


where 


(60) 


*  Re[4fn]  :  UP  =  Im[A#„] 

Therefore,  with  all  terms  to  order  A«f  retained  in  iln,  and  assuming  A^°^0. 

E&*Ht ;«*)=  -  g— , ' «/  [t [( 1  +a,)+d2(f 2 


where 


A*«  [Ac**]  1  ^2)(Ao|] 

5l  =  vij0>  1  2/1  j  j1  +  4  >^(0)  2/i 

-ism 


It  is  evident  from  (81)  that  the  envelope  functions  are,  nearly,  simple  Gaus¬ 
sian  exponentials  with  a  constant  coefficient  equal  to  the  spectral  amplitude  of 
the  signal  atu  =  ut.  The  departure  from  this  simple  form  involves  modulation  by 
the  square  root  of  a  polynomial  function  of  ( t  where  in  the  case  when 

is  truncated  after  terms  of  order  A&>£/2/?,  this  polynomial  is  of  second  order. 
For  pulse-like  signals  conforming  with  our  criteria  of  smoothness,  that  is  with 

|4i°W|  >  |*!”(  uk  )|>l4*'(«.)|>- 

and  for  narrow-band  filters  designed  according  to  the  criteria  of  (28)  or  (29), 


then 


\6g\  <  |dj|  «  1 


and  the  effect  of  this  extra  modulation  term  is  small. 

The  instantaneous  phase  for  a  single  pulse  is,  from  the  definitions  in  (59) 


and  the  relations  (57)  and  (58), 


(63) 


xln)  =  tan*1 
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t>k 


is_ 

2  ok 


with  and  tf  defined  previously.  This  result  is.  of  course,  accurate  when 

t  is  close  to  tjn\  the  group  arrival  time  for  the  nfh.  signal  pulse.  When  t  is,  in 
fact,  equal  to  tjn\  then  the  phase  factor  is  negligible,  as  can  be  seen  from 
(56).  Since  xLn^  Is  8  calculable  function  of  time,  then  (63)  can  be  used  to  com¬ 
pute  the  phase  time  t£n);  given  that  the  group  time  can  be  otherwise  deter¬ 
mined  (see  the  next  section)  and  the  phase  factor  is  either  entirely  negligi¬ 
ble,  which  is  usually  the  case,  or  has  been  estimated  from  the  variation  of  the 
measured  t versus  frequency.  In  this  case,  f^n)  can  be  used  to  obtain  a 
phase  velocity  estimate,  as  a  function  of  frequency,  using  its  definition  given  in 
(15).  In  addition,  the  instantaneous  phase,  when  measured  for  three  com¬ 
ponents  of  the  displacement,  can  be  used  to  determine  the  polarization  of  the 
wave  field  as  a  function  of  frequency.  In  all  cases,  the  functions  derived  are 
valid  approximations  to  the  true  variable,  such  as  phase  velocity  and  polariza¬ 
tion.  only  at  the  times  very  near  the  group  times,  t§ ,  for  the  signal  in  question. 

The  instantaneous  frequency  for  the  nth  signal  pulse  flin\  is,  from  (51): 

0in,(<  :ok)  *  =  [&  ~  5.  ^r]/  [5?  +  S?]  (64) 


The  time  derivatives  appearing  here  are,  using  (57)  and  (58): 
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Forming  the  relation  for  f)J*)  therefore  gives,  from  (64): 
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where  the  terms  involving  the  time  factor  (f  -  </n^),  appearing  in  the  relations 
of  (65),  cancel  out. 

From  (66)  it  is  apparent  that,  to  first  order  the  instantaneous  frequency 
function  is  constant  and  equal  to  the  center  frequency  of  the  Gaussian  filter 
used.  However,  accounting  for  terms  of  second  order  in  the  filter  bandwidth 
Aw*.  produces  the  interesting  result  that  the  constant  frequency  value  is 
shifted  somewhat  and  in  a  manner  proportional  to  the  first  frequency  derivative 
of  the  amplitude  spectrum  of  the  signal.  More  specifically,  when  the  ratio  of 
the  first  frequency  derivative  of  the  spectrum  to  the  spectral  amplitude  itself 
(Le.,  ylr^,VA^0,)  becomes  relatively  large,  as  would  occur  in  the  vicinity  of  a 
sharp  spectral  minimum,  then  this  term  can  become  quite  significant  in  size 
and  result  in  a  large  deviation  from  the  filter  center  frequency,  uk.  In  particu¬ 
lar  this  second  order  term  should  give  a  double  spike  similar  to  a  derivative  of  a 
delta  function  centered  at  the  spectral  minimum  of  the  signal.  Such  behavior 
is,  in  fact,  observed  to  occur,  as  is  noted  in  the  examples  of  later  sections. 

Relationships  between  Envelope  Function  Maxima  and 
Signal  Group  Arrival  Times 

The  most  important  and  useful  feature  of  the  envelope  function  is  the  fact 
that  its  maxima  occur  at,  or  more  accurately  near,  the  times  t  =t}n\  This  allows 
the  group  arrival  times  to  be  estimated  for  each  frequency  ok.  Further  at  these 
times  the  amplitude  of  the  envelope  will  be  near 

To  show  this  explicitly  and  to  obtain  quantitative  results,  consider  the  loca¬ 
tion  of  the  envelope  function  maxima  at  t  as  given  by 


We  have,  using  (61)  and  (62): 
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Thus  maxima  occur  at: 
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with  amplitudes  given  by 
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when  is  represented  by  (58).  That  is,  when  terms  of  order  [Auf/20]z  are 
negligible  in  the  modulation  function  JU^,  as  in  (54),  then  envelope  maxima  are 
predicted  to  occur  at  precisely  the  t ^  group  times.  Further  the  envelope 
amplitude,  at  its  maxima,  gives  a  measure  of  the  signal  Fourier  spectrum  at 
each  frequency  &>*.  Note  that  the  maxima  times  do  not  explicitly  depend  on  the 
condition  that  polynomial  coefficients  Sx  and  d2  be  small  -  although  for  a  trun¬ 
cated  polynomial  to  be  a  valid  approximation  for  Mn  it  is  necessary  that  this  be 
so.  Clearly  at  some  stage  the  "roughness"  in  the  amplitude  and  group  arrival 
time  as  functions  of  frequency  must  result  in  a  shift  of  the  maxima  away  from 
the  t values. 

In  fact,  it  is  not  difficult  to  obtain  estimates  of  the  shift  in  the  envelope 
maxima  due  to  non-negligible  effects  of  dispersion  and  amplitude  spectra 
roughness  (  as  measured  by  the  frequency  derivatives  of  the  signal  amplitude 
spectra  and  group  arrival  times).  This  is  simply  achieved  by  including  the  next 
higher  order  terms  for  the  representation  of  and  in  the  expression 
(57),  for  the  envelope  function.  In  particular,  we  have: 
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Now  forming  the  derivative  dEir*/  dt  and  retaining  terms  to  order  but.  we  find 
maxima  of  E&n^  at  times  f  given  by  the  zeros  of  : 


tfo  +  «i(‘4r)-*,(n,)  =  o 
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With  terms  of  order  but  retained,  these  coefficients  are: 
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Thus  the  maxima  of  the  envelope  function  are  at  times  fm  given  by: 


U")  =  +  a<£*> 


where 


sor][-57-J,.t3|-^-i,.|ir/ 11 + HsH'lsH  W. 


It  is  sufficient,  for  "well  designed"  narrow-band  Gaussian  filters  operating  on 
pulse-like  signals,  to  retain  only  terms  of  order  Awf  in  so  that  in  this  case. 
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+  O(Acjjf) 


The  form  of  6t^  is  such  that  the  shifts  in  the  envelope  function  maxima  away 
from  f^m)  values  are  small  for  pulse-like  signals,  for  which: 
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if  the  filter  design  criteria  of  (28)  or  (29)  are  satisfied.  If  the  signals  in  question 
do  not  conform  to  the  spectral  smoothness  criteria  above,  then  other  design 
criteria  based  on  the  reduction  of  ttjjj?*  to  a  small  number  (  relative  to  tjn)), 
would  be  required.  In  either  case,  equation  (72)  or  (73)  can  be  used  to  (itera¬ 
tively)  correct  the  estimates  of  using  measured  tj?)  values  and  first  order 
estimates  of  the  various  derivatives  appearing  in  these  equations. 
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Since  it  is  *4,*)  that  is  measurable,  then  clearly  it  is  the  amplitude  of  the 
envelope  function  at  l=<4P*  that  can  provide  the  consistently  measured  esti¬ 
mate  of  the  spectra  amplitude  of  the  signal.  This  peak  value  for  at  t  is, 
from  (69)  and  subsequent  relations. 


1  +  & 

Aw|| 

20 

(74) 


when  carried  to  terms  of  order  Awf.  Comparison  with  the  results  for  £jfn)  at 
t  ~t^n\  given  by  equation  (65),  shows  that  we  get  the  same  value  for  the 
envelope  maxima  at  <4P*  when  terms  of  order  Aw*  are  neglibible.  Thus  the  small 
shift  in  time  of  the  envelope  peak  from  the  group  arrival  time  of  the  signal  does 
not  significantly  change  the  envelope  amplitude  value,  since  differences  in  the 
envelope  amplitudes  at  (#W  and  +  are  of  order  Aw*. 

The  value  of  measured  at  will  be  very  nearly  equal  to  A^\uk)  - 
the  Fourier  spectral  amplitude  of  the  signal  at  w*  -  when  the  design  criteria  of 
(28)  or  (29)  are  met  -  since  the  second  term  in  (74)  will  be  much  less  than 
unity.  We  can  therefore  estimate  the  signal  spectrum  using 


4P»(w*) 


as  a  good  first  approximation,  and  then  obtain  First  and  second  derivative  esti¬ 
mates  from  differencing  the  41°)  values  obtained  at  the  set  of  filter  center  fre¬ 
quencies  (w*).  These  derivative  estimates  can  then  be  used  to  obtain  a  higher 
order  iterative  estimate  of  A$?\  employing  (74).  This  procedure  can,  of  course, 
also  be  used  to  get  higher  order  estimates  of  the  values,  using  (73)  and  (74) 
together. 
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CHAPTER  I 


INTRODUCTION 

Historically,  normal  mode  theory  has  been  restricted  in  its  appli¬ 
cations  to  low  frequency  bandwidths;  namely  to  low  frequency  spherical 
earth  normal  modes  and  to  Rayleigh  and  Love  surface  waves  for  flat  layered 
structures.  In  this  dissertation  I  will  show  how  spectral  solutions  of  the  elas¬ 
tic  wave  equation  can  be  used  to  compute  complete,  high  frequency  syn¬ 
thetic  seismograms  for  flat,  plane  layered,  and  laterally  homogeneous  struc¬ 
tural  models  in  an  efflcient  manner.  The  method  which  I  developed  is  most 
useful  for  computing  synthetic  seismograms  in  the  zero  to  ten  Hertz  fre¬ 
quency  range  and  for  source-receiver  distances  of  10  to  1000  km.  I  have  also 
been  able  to  successfully  apply  the  method  to  exploration  geophysics  prob¬ 
lems  with  frequency  bandwidths  of  100  Hertz  and  source-receiver  distances 
of  several  km.  By  making  a  simple  modification  to  the  structural  model  I 
am  able  to  use  this  method  to  compute  the  P  and  S  body  waves,  in  addition 
to  the  surface  waves,  using  only  normal  modes  with  real  eigenwavenumbers. 
I  am  thus  able  to  approximate  a  complete  solution  of  the  elastic  wave  equa¬ 
tion  with  a  mode  sum  which  makes  this  approach  much  more  efflcient  than 
existing  complete  solution  methods  which  are  based  upon  direct  numerical 
integration  such  as  the  reflectivity  method. 

A  variety  of  approachs  have  been  used  to  synthesize  the  P  and  S 
body  waves  for  laterally  homogeneous  structures,  and  all  of  these  approachs 
start  with  a  doubly  transformed  version  of  the  elastic  wave  equation  which 
remove  the  derivatives  of  time  and  the  horizontal  space  coordinates.  The 
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methods  differ  by  the  solution  of  the  transformed  wave  equation  which  is 
used  and  by  the  way  in  which  the  two  inverse  transform  integrals  are 
evaluated.  The  asymptotic  methods  evaluate  either  one  or  both  of  these 
integrals  analytically  by  approximating  the  integrand  function  with  suitable 
asymptotic  expansions,  and  the  resulting  solution  is  a  decomposition  of  the 
complete  solution  in  terms  of  rays.  These  ray  theories  work  well  for  syn¬ 
thesizing  particular  phases,  but  they  can  be  cumbersome  to  use  when  trying 
to  compute  the  complete  solution,  especially  in  certain  distance  ranges.  For 
very  near  field  distances  (zero  to  ten  km)  and  at  teleseismic  distances 
(greater  than  1000  km)  and  for  sources  at  typical  earthquake  depths,  the 
elastic  energy  is  propagating  in  a  fundamentally  vertical  direction  in  the 
crust.  This  means  that  the  scattering  off  of  the  large  discontinuities  in  the  P 
and  S  wave  velocities  that  can  occur  in  the  crust  and  at  the  Mohorovicic 
discontinuity  (Moho)  can  be  represented  by  a  relatively  small  number  of  ray 
paths,  and  so  the  major  applications  of  ray  theory  have  been  in  these  dis¬ 
tance  ranges.  In  the  10  to  1000  km  distance  range  however  the  crust  acts 
as  a  waveguide  and  most  of  the  elastic  energy  is  contained  within  this 
waveguide  and  propagates  in  a  fundamentally  horizontal  direction.  In  this 
distance  range, 'when  using  a  detailed  crust  and  upper  mantle  model,  ray 
theories  require  the  a  priori  specification  of  a  very  large  number  of  ray 
paths  to  synthesize  the  complete  solution  for  the  P  and  S  body  waves  (a 
typical  example  of  this  is  the  Pg  coda  which  is  seen  in  the  western  United 
States). 

Another  category  of  seismogram  synthesis  techniques  which  have 
been  used  over  the  past  fifteen  years  are  the  complete  solution  methods  and 
these  methods  all  differ  from  the  ray  theoretical  methods  by  the  solution  of 
the  doubly  transformed  elastic  wave  equation  which  is  used.  The  complete 
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solution  methods,  as  the  name  implies,  use  a  complete  solution  of  the  dou¬ 
bly  transformed  elastic  wave  equation  by  allowing  for  all  possible  P  and  S 
wave  propagation  throughout  the  structural  model  and  by  treating  this  as  a 
boundary  value  problem.  When  one  uses  a  complete  solution  method  it  is 
only  necessary  to  specify  the  structural  model  and  frequency  bandwidth  to 
compute  synthetic  seismograms  and  these  synthetics  contain  an  infinity  of 
rays.  In  contrast,  ray  theoretical  methods  require  the  user  to  specify  ray 
paths  and,  based  upon  these  specified  ray  paths,  an  incomplete  solution  is 
obtained.  It  is  the  use  of  this  incomplete  solution,  in  addition  to  certain 
other  approximations  which  are  usually  made,  that  causes  the  ray  theoreti¬ 
cal  methods  to  be  much  more  efficient  than  complete  solution  methods  and 
it  is  primarily  this  efficiency  that  has  made  ray  theoretical  methods  so 
popular. 

The  complete  solution  methods  are  themselves  broken  down  into 
two  general  categories  which  I  refer  to  as  the  reflectivity  method  and  the 
spectral  method  and  these  methods  differ  in  the  way  in  which  the  two 
inverse  transform  integrals  are  evaluated.  I  am  using  the  reflectivity  method 
to  refer  to  all  methods  which  compute  both  integrals  in  a  direct  numerical 
manner  although  the  original  reflectivity  method,  as  developed  by  Fuchs 
and  Muller  (1971),  is  rather  restrictive  in  terms  of  the  horizontal  phase 
velocity  range  over  which  it  works.  Owing  to  recent  developments  (Kind 
(1978),  Kennett  and  Kerry  (1979),  Kennett  (1980),  Cormier  (1980), 
Bouchon  (1981))  the  reflectivity  method  can  now  be  used  to  compute  com¬ 
plete  seismic  codas  for  arbitrary  frequency  bandwidths  and  source  depths 
and  for  vertically  inhomogeneous  structural  models.  If  the  reflectivity 
method  were  also  efficient,  then  the  problem  of  computing  complete  syn¬ 
thetic  seismograms  for  flat,  laterally  homogeneous  earth  models  could  be 
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considered  to  be  solved  however  the  reflectivity  methods  are  very  computa¬ 
tionally  expensive  and  this  expense  is  proportional  to  the  product  of  the  fre¬ 
quency  bandwidth  and  the  maximum  source-receiver  distance.  In  fact  the 
computational  expense  is  so  high  that  it  limits  the  use  of  these  methods  to 
low  frequency  bandwidths  for  cases  where  the  complete  seismic  coda  is  com¬ 
puted. 

One  obvious  way  of  increasing  the  efficiency  of  complete  solution 
methods  is  to  analytically  evaluate  at  least  one  of  the  inverse  transform 
integrals  as  is  done  by  most  ray  theoretical  methods.  Unfortunately  the 
complexity  of  the  complete  solution  integrand  function,  for  generally  com¬ 
plex  structural  models,  frustrates  efforts  to  apply  the  types  of  techniques 
which  are  used  by  ray  theoretical  methods  to  eliminate  the  numerical 
evaluation  of  the  inverse  transform  integrals.  There  is  one  straightforward 
method  however  which  we  can  always  use  to  evaluate  at  least  part  of  one  of 
the  inverse  transform  integrals  analytically  and  this  method  makes  use  of 
the  residue  theorem.  We  can  extend  the  integration  path  in  the  complex 
plane  into  a  closed  contour  and  then  evaluate  the  original  integral  in  terms 
of  a  sum  of  residues  which  are  caused  by  the  poles  of  the  integrand  function 
which  are  contained  within  the  contour  of  integration.  These  poles 
correspond  to  flat  earth  normal  modes  and  for  structural  models  which  have 
totally  reflective  top  and  bottom  boundaries,  one  of  the  inverse  integral 
transforms  can  be  expressed  as  an  infinite  sum  of  normal  modes.  The  cases 
of  most  interest  in  seismology  are  for  structural  models  which  have  a  free 
top  boundary  and  an  infinite  homogeneous  half  space  on  the  bottom  and  for 
these  situations  one  of  the  inverse  integral  transforms  can  be  expressed  as  a 
finite  sum  of  normal  modes  along  with  branch  cut  integrals  which  come 
about  due  to  the  semi-infinite  nature  of  the  structural  model. 
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1  refer  collectively  to  methods  which  use  the  residue  theorem  to 
evaluate  one  of  the  inverse  transform  integrals  as  the  spectral  method. 
Although  this  method  has  not  been  as  popular  as  the  reflectivity  method  for 
computing  complete  synthetic  seismograms,  it  has  undergone  parallel 
developments  and  improvements  with  the  reflectivity  method.  The  first  uses 
of  the  spectral  method  actually  predate  the  reflectivity  method  since  the 
spectral  method  is  the  basis  behind  the  computation  of  synthetic,  flat  earth, 
Rayleigh  and  Love  fundamental  surface  waves  however,  among  the  earliest 
uses  of  the  spectral  method  for  computing  a  substantial  portion  of  the  com¬ 
plete  seismic  coda,  are  those  of  Knopoff,  et.  al.  (1973),  Nakanishi,  et.  al. 
(1977),  and  Mantovani,  et.  al.  (1977)  who  used  a  sum  of  SH  normal  modes 
to  compute  synthetic  SH  seismograms.  More  recently  Swanger  and  Boore 
(1978)  computed  both  SH  and  P-SV  synthetic  seismograms  for  near  field 
strong  motion  studies  using  a  normal  mode  sum.  All  of  these  uses  of  the 
spectral  method  had  one  thing  in  common  which  was  that  a  small  number 
of  normal  modes  was  included  in  the  mode  sum  and  this  resulted  in  rather 
incomplete  solutions  to  the  elastic  wave  equation. 

One  fundamental  difficulty,  which  had  a  large  effect  on  both  the 
development  of  the  complete  solution  methods  and  their  ranges  of  applica¬ 
bility,  was  a  numerical  instability  which  seened  to  be  inherent  in  the  com¬ 
plete  solution  form  of  the  P-SV  integrand  function.  This  numerical  instabil¬ 
ity  was  always  associated  with  the  presence  of  inhomogeneous  or  evanescent 
waves  within  the  elastic  medium  which  exist  at  horizontal  phase  velocities 
that  are  less  than  the  local  P  or  S  wave  velocity.  This  problem  was  first 
recognized  by  Dorman,  et.  al.  (1960)  when  the  complete  matrix  solution  for¬ 
malism  of  Thomson  (1950)  and  Haskell  (1953)  was  applied  to  the  problem 
of  computing  Rayleigh  dispersion  curves  using  an  early  digital  computer. 
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Although  the  numerical  problem  was  circumvented  for  the  computation  of 
the  fundamental  Rayleigh  dispersion  curve  by  the  use  of  structural  layer 
reduction,  it  wasn't  until  Knopoff  (1964)  applied  Laplace's  development  by 
minors  to  the  problem  of  computing  the  Rayleigh  characteristic  function 
that  the  numerical  problem  was  formally  addressed  in  an  analytical  manner. 
Knopoffs  work  was  the  seed  for  an  area  of  research  which  was  followed  by 
Dunkin  (1965),  Watson  (1970)  and  most  recently  by  Abo-Zena  (1979)  and 
all  of  this  research  was  aimed  at  streamlining  KnopoiTs  original  method  and 
understanding  why  it  worked  as  well  as  it  did.  The  method  completely 
solved  the  numerical  instability  problem,  as  related  to  computing  the  Ray¬ 
leigh  characteristic  function  and  thus  the  Rayleigh  dispersion  curves,  for  all 
frequency  bandwidths  and  structural  models. 

Although  KnopofTs  method  works  quite  well  for  computing  the 
Rayleigh  characteristic  function  it  does  not  address  the  problem  of  numeri¬ 
cal  instabilities  which  are  present  in  the  computation  of  the  depth  depen¬ 
dent  stress  and  displacement  eigenfunctions.  Because  of  this  it  has  not  been 
possible  to  use  the  spectral  method  to  compute  P-SV  synthetic  seismograms 
for  buried  sources  and  at  high  frequencies  when  using  the  Knopoff  modified 
version  of  the  Thomson-Haskeil  matrix  formalism.  Thus  we  can  see  that 
there  are  actually  two  numerical  instability  problems  which  must  be  solved 
in  order  to  use  the  spectral  method  to  compute  P-SV  synthetic  seismograms 
for  arbitrary  sources  and  frequencies,  and  my  research  has  focused  on 
extending  Knopoffs  method  to  solve  the  eigenfunction  numerical  problem. 

The  same  numerical  difficulties  that  plagued  the  spectral  method 
were  encountered  with  the  reflectivity  method.  The  original  method  as 
given  by  Fuchs  and  Muller  was  based  directly  on  the  original  Thomson- 
Haskeil  matrix  method  and  they  avoided  the  numerical  problems  by 
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limiting  the  range  of  phase  velocities  to  those  for  which  there  were  no 
evanescent  waves  within  the  structure.  This  narrowed  the  range  of  applica¬ 
bility  of  the  reflectivity  method  until  Kind  (1978)  applied  KnopofTs  method 
to  remove  the  numerical  instabilities.  Kind  also  reformulated  the  integrand 
function  in  a  manner  which  he  claimed  eliminated  the  numerical  instabili¬ 
ties  associated  with  buried  sources.  A  completely  different  and  novel 
approach  was  taken  by  Kennett,  et.  al.  (1978),  Kennett  and  Kerry  (1979) 
and  Kennett  (1980)  to  solve  the  numerical  problems  of  the  reflectivity 
method  which  was  not  an  extension  of  KnopofTs  original  work.  They 
showed  how  the  doubly  transformed  complete  solution  of  the  elastic  wave 
equation  could  be  expressed  in  terms  of  a  set  of  generalized  reflection  and 
transmission  functions.  They  then  showed  how  all  growing  exponential 
solutions  could  be  eliminated  in  this  ray-like  representation  of  the  complete 
solution  and  this  eliminated  the  numerical  instabilities  for  all  source  depths. 

The  most  recent  developments  of  the  spectral  method  have 
focused  on  solving  the  numerical  problems  for  arbitrary  frequencies,  source 
depths,  and  structural  models  and  extending  the  completeness  of  the  spec¬ 
tral  solutions.  All  of  the  earlier  uses  of  the  spectral  method  used  only  a 
small  number  of  normal  modes  and  thus  produced  rather  incomplete  solu¬ 
tions.  Harvey  (1981)  and  Kerry  (1981)  were  the  first  researchers  to  use  all 
of  the  normal  modes  w'ith  real  eigenwavenumbers  and  I  refer  to  spectral 
synthetic  seismograms  produced  in  this  manner  as  locked  mode  synthetics. 
The  branch  cut  integral  contributions  are  ignored  when  using  the  locked 
mode  method  which  gives  the  most  complete  spectral  solution  possible 
without  using  numerical  integration  or  without  locating  poles  off  of  the  real 
wavenumber  axis.  At  about  the  same  time  Wang  and  Herrmann  (1980) 
developed  a  truly  complete  spectral  solution  by  including  both  a  numerical 
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integration  of  the  branch  cuts  along  with  all  of  the  locked  mode  residues. 

As  with  the  reflectivity  method,  the  development  of  numerically 
stable  algorithms  to  use  in  the  spectral  method  has  followed  two  paths. 
Kerry  (1981)  directly  applied  the  methods  developed  by  Kennett,  et.  al. 
which  reformulated  the  frequency-wavenumber  solution  of  the  elastic  wave 
equation  in  terms  of  generalized  reflection  and  transmission  functions.  I 
chose  to  extend  Knopoffs  method,  which  operates  directly  with  the 
Thomson-Haskell  matrix  formalism  and  uses  stress  and  displacement  func¬ 
tions  to  describe  the  elastic  propagation,  and  the  numerical  algorithms 
which  I  developed  constitute  a  substantial  portion  of  my  research.  Kerry's 
locked  mode  method  comes  the  closest  to  my  own  both  in  terms  of  the  basic 
way  in  which  it  works  and  its  range  of  applicability  and  1  will  be  comparing 
these  two  methods  throughout  this  dissertation. 

The  research  which  I  will  be  describing  in  the  following  chapters 
has  as  its  basic  objective  the  development  of  an  optimally  fast,  accurate  and 
complete  spectral  method  for  computing  P-SV  synthetic  seismograms  for 
flat,  plane  layered,  laterally  homogeneous  structural  models.  An  equally 
important  criterion  which  I  placed  on  the  method  is  that  it  work  for  the 
widest  possible  range  of  frequency  bandwidths  and  structural  models.  I 
adopted  this  last  criterion  to  cover  problems  such  as  high  frequency  P,  and 
SB  wave  propagation  in  oceanic  structural  models  with  liquid  and  near¬ 
liquid  layers  and  very  high  frequency  near  surface  wave  propagation  in 
structural  models  which  have  numerous  strong  low  velocity  zones  for  prob¬ 
lems  in  earthquake  hazards  engineering  and  geophysical  exploration. 

Since  eigensolutions  of  the  elastic  wave  equation  can  be  thought 
of  as  providing  optimal  sampling  in  the  wavenumber  domain,  spectral 
methods  should  be  able  to  provide  the  most  efficient  way  to  compute 


complete  synthetic  seismograms.  1  feel  that  one  of  the  fundamental  reasons 
that  the  spectral  has  not  seen  the  popularity  of  the  reflectivity  method  has 
been  the  lack  of  an  efficient  and  reliable  normal  mode  searching  algorithm. 
A  substantial  portion  of  my  research  has  been  devoted  to  developing  a  truly 
fast,  accurate  and  completely  reliable  mode  searching  algorithm  in  order  to 
realize  the  potential  efficiency  of  the  spectral  method.  Also,  in  order  to 
maximize  the  efficiency,  it  was  desirable  to  compute  nearly  complete  syn¬ 
thetic  seismograms  without  making  numerical  evaluations  of  the  branch  cut 
integrals.  For  normal  structural  models  locked  mode  synthetic  seismograms 
will  not  contain  any  P  wave  arrivals  since  they  are  typically  part  of  the 
branch  cut  integral  contribution  however,  by  making  a  simple  modification 
to  the  structural  model,  one  can  significantly  extend  the  phase  velocity 
range  which  will  be  represented  in  the  locked  mode  synthetic  seismograms 
while  simultaneously  insuring  that  a  certain  time  window  within  the  syn¬ 
thetic  seismograms  will  be  uncontaminated  by  the  structural  modification. 
When  the  locked  method  is  applied  to  such  modified  structural  models,  I 
refer  to  this  as  the  locked  mode  approximation  which  produces  nearly  com¬ 
plete  synthetic  seismograms  while  only  using  normal  modes  with  real 
eigenwavenumbers . 

In  chapter  two  I  review  the  basic  theory  for  computing  flat  earth 
normal  modes  and  the  resulting  displacements.  Most  of  this  draws  upon 
previously  published  work  however  1  will  present  a  complete  and  consistent 
derivation  starting  with  the  elastic  wave  equation  and  ending  with  the  spec¬ 
tral  solution  for  flat,  plane  layered,  isotropic  and  laterally  homogeneous 
models.  During  this  derivation  I  will  indicate  the  departure  points  of  the 
various  seismogram  synthesis  methods  and  I  will  also  point  out  the  sources 
of  the  numerical  instabilities.  The  final  solution  which  I  derive  will  be 
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expressed  in  terms  of  proper  and  improper  eigenfunctions  representing 
respectively  the  discrete  normal  modes  and  the  continuous  branch  cut 
integrals. 

In  chapter  three  I  will  address  the  numerical  instabilities  which 
are  inherent  in  the  Thomson-Haskell  matrix  formulation  and  1  will  present 
two  methods  which  overcome  these  instabilities.  I  start  by  reviewing 
Knopoff  s  method,  along  with  the  later  work  which  was  based  upon  it,  and  I 
give  the  derivation  of  a  numerically  stable  computation  of  the  Rayleigh 
characteristic  function.  I  then  describe  a  simple  method  for  stabilizing  the 
depth  dependent  stress  and  displacement  eigenfunctions  which  works  quite 
well  for  completely  elastic  structural  models  with  monotonic  ally  increasing 
velocities  with  depth  except  for  no  more  than  one  weak  low  velocity  zone. 
The  locked  mode  method  when  using  this  algorithm  is  functionally 
equivalent  to  the  method  developed  by  Kerry  and  I  discuss  the  restrictions 
which  these  methods  have  in  common.  I  then  proceed  to  derive  a  much 
more  robust  method  of  computing  numerically  stable  eigenfunctions  which 
has  virtually  no  restrictions.  This  method  works  for  oceanic  models  as  well 
as  for  complex  elastic  models  with  multiple  strong  low  velocity  zones  and  at 
arbitrary  frequency  b&ndwidths.  I  finally  show  numerical  examples  of  depth 
dependent  eigenfunctions  using  the  two  methods  for  several  different  struc¬ 
tural  models. 

In  chapter  four  I  describe  the  numerical  algorithms  and  computer 
programs  which  I  developed  to  implement  the  computation  of  locked  mode 
synthetic  seismograms  from  a  starting  structural  model  to  the  final,  three 
component  time  traces  at  specified  receiver  locations.  The  first  step  in  this 
process  is  the  normal  mode  searching  algorithm  which  locates  the  Rayleigh 
and  Love  dispersion  curves  and  I  go  into  considerable  detail  to  describe  the 


algorithm  which  I  developed.  The  next  step  involves  the  computation  of 
certain  partial  derivatives  which  are  necessary  for  eigenfunction  normaliza¬ 
tion  and  which  I  use  in  a  first  order  perturbation  approximation  to  account 
for  the  effects  of  frequency  dependent  anelastic  attenuation  in  the  structural 
model.  I  give  analytic  expressions  for  these  derivatives  and  I  show  how  the 
first  order  attenuation  approximation  can  be  computed  and  applied  in  a  fast 
manner  which  does  not  require  the  use  of  complex  arithmetic.  In  the  next 
section  I  describe  the  actual  computer  programs  which  I  wrote  and  the  user 
interface  with  these  programs.  I  discuss  practical  matters  such  as  the  algo¬ 
rithms  which  require  the  use  of  double  precision  arithmetic,  the  amount  of 
core  memory  required  by  each  program,  the  approximate  run  times  of  the 
programs,  and  the  structure  and  size  of  the  intermediate  data  files  which 
link  the  programs  together.  I  then  formally  present  the  locked  mode  approx¬ 
imation  and  show  when  the  approximation  breaks  down.  In  this  section  I 
also  show  synthetic  seismograms  produced  by  the  locked  mode  approxima¬ 
tion  and  how  spurious  arrivals  caused  by  the  approximation  can  be  con¬ 
trolled.  Finally  I  show  comparisons  of  synthetic  seismograms  produced  by 
the  locked  mode  approximation  with  synthetics  for  the  same  structural 
models  which  were  generated  using  other  complete  solution  methods. 

Chapter  five  is  devoted  to  showing  examples  of  synthetic  seismo¬ 
grams  produced  by  the  locked  mode  approximation  for  a  variety  of  fre¬ 
quency  bandwidths,  structural  models  and  source-receiver  distances.  In  the 
first  part  of  the  chapter  I  show  a  number  of  examples  which  illustrate  the 
characteristics  of  elastic  wave  propagation  which  can  be  seen  when  using  a 
complete  solution  method.  I  then  show  an  example  in  which  theoretical 
seismograms  using  the  locked  mode  approximation  are  compared  to  real 
observed  data.  This  example  involves  modelling  an  underground  nuclear 


explosion  which  took  place  in  northern  New  Mexico  and  for  which  receivers 
recorded  seismograms  in  the  10  to  100  km  distance  range  and  the  0  to  10 
Hertz  frequency  range.  This  is  a  well  constrained  example  since  both  the 
source  and  the  structure  were  known  and  since  the  structure  was  closely 
approximated  by  a  flat  layered  half  space. 

In  chapter  six  I  conclude  the  dissertation  by  summarizing  the 
relative  advantages  and  disadvantages  of  the  locked  mode  approximation 
when  compared  to  other  seismogram  synthesis  techniques.  The  possibility  of 
future  extensions  of  my  research  are  discussed  which  would  extend  the 
range  of  applicability  of  the  methods  which  I  have  developed  while  main¬ 
taining  the  efficiency. 


CHAPTER  H 

THE  FUNDAMENTALS  OF  NORMAL  MODE  THEORY 
FOR  FLAT  LAYERED  STRUCTURES 

Although  the  development  of  the  theoretical  basis  for  elastic  wave 
propagation  can  be  traced  all  the  way  back  to  Lord  Rayleigh’s  time,  the 
first  occurrence  of  the  complete  solution  of  the  elastic  wave  equation  in 
cylindrical  coordinates  and  for  flat,  plane  layered,  isotropic  and  laterally 
homogeneous  elastic  media  is  given  by  Sezawa  (1931).  He  uses  a 
transformed  version  of  the  elastic  wave  equation  which  eliminates  all  deriva¬ 
tives  except  for  the  depth  derivatives  and  thus  he  reduces  the  problem  of 
solving  the  elastic  wave  equation  to  one  of  solving  several  ordinary  differen¬ 
tial  equations  and  evaluating  inverse  integral  transforms.  The  form  of  the 
solution  which  Sezawa  presents  is  identical  to  that  used  by  all  modern  day 
seismologists  who  work  in  cylindrical  coordinates. 

Sezawa  did  not  give  solutions  for  the  depth  dependent  ordinary 
differential  equations  for  arbitrary  vertical  velocity  distributions  and  the 
next  major  development  was  aimed  at  solving  these  equations.  Thomson 
(1950)  derived  a  solution  for  the  depth  dependent  differential  equations 
which  was  expressed  in  terms  of  recursive  matrix  multiplications. 
Thomson’s  method  applied  to  elastic  structures  which  were  composed  of 
planar,  isotropic  and  completely  homogeneous  layers  which  were  welded 
together  along  their  top  and  bottom  surfaces.  Although  each  layer  had  to  be 
homogeneous  there  was  no  restriction  on  the  number  of  different  layers  that 
were  welded  together  or  on  the  thinness  of  each  layer  and  so,  in  the  limit  of 
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an  infinite  number  of  infinitesimally  thin  layers,  Thomson's  method  pro¬ 
vided  complete  solutions  of  the  depth  dependent  differential  equations  for 
arbitrary  vertical  velocity  profiles.  Haskell  (1953)  embraced  Thomson’s 
approach  and  basically  streamlined  it  and  applied  it  to  produce  numerical 
computations  of  the  fundamental  Rayleigh  dispersion  curve  for  several  sim¬ 
ple  structural  models.  The  Thomson- Haskell  matrix  method  has  provided  a 
formalism  for  obtaining  complete  solutions  of  the  elastic  wave  equation  and, 
as  with  Sezawa’s  work,  this  has  become  one  of  the  basic  tenets  of  a  large 
branch  of  subsequent  research. 

It  wasn’t  until  digital  computers  became  generally  available  that 
the  next  major  developments  in  the  theory  of  seismic  wave  propagation  in 
cylindrical  coordinates  took  place.  In  the  early  sixties  a  number  of  research¬ 
ers  applied  the  Thomson-Haskeil  matrix  method  to  write  computer  pro¬ 
grams  for  locating  Rayleigh  dispersion  curves  (e.g.  Dorman,  et.  al.  (1960), 
Press,  et.  al.  (1961)).  Numerical  instabilities  in  the  Thomson-Haskeil 
method  were  discovered  at  this  time  however  these  instabilities  occurred  at 
frequencies  which  were  above  the  range  of  interest  of  the  researchers.  Har- 
krider  (1964)  made  the  next  major  contribution  by  showing  how  to  compute 
time  domain  synthetic  seismograms  for  Rayleigh  and  Love  surface  waves 
and  for  arbitrary  point  sources  at  arbitrary  depths.  This  research  was  based 
directly  on  the  Thomson-Haskeil  matrix  method  and  it  suffered  an  addi¬ 
tional  numerical  instability  associated  with  the  source  depth. 

In  the  following  chapter  I  give  a  condensed  yet  complete  sum¬ 
mary  of  linear  elastic  wave  theory  for  flat,  plane  layered,  laterally  homo¬ 
geneous  structures.  This  treatment  will  draw  primarily  from  the  work  of 
Harkrider  (1964)  and  Ben-Menahem  and  Singh  (1972  and  1981)  who  show 
how  a  complete  solution  of  the  elastic  wave  equation  can  be  expressed  in 


terms  of  vector  cylindrical  harmonics,  but  in  order  to  preserve  continuity,  I 
will  use  my  own  notation  and  layering  conventions.  I  will  show  how  normal 
mode  solutions  are  derived  from  the  general  solutions  of  the  elastic  wave 
equation  and  to  what  extent  these  solutions  are  related  to  other  types  of 
solutions  such  as  ray  theoretical  and  reflectivity  solutions.  I  will  also  point 
out  both  accuracy  and  efficiency  problems  which  come  up  in  the  numerical 
implementation  of  the  theory.  These  problems  will  be  addressed  in  detail  in 
subsequent  chapters  and  constitute  the  major  portion  of  my  research. 

2.1  Basic  Assumptions 

In  order  to  obtain  tractable  solutions  of  the  elastic  wave  equation 
it  is  necessary  to  make  a  number  of  assumptions.  Most  of  the  assumptions 
I  have  made  fall  in  this  category  and  are  absolutely  necessary  in  order  to 
solve  the  problem  at  all.  I  have  however  made  a  few  additional  assump¬ 
tions  which  were  made  primarily  to  narrow  the  range  of  the  problem  and  I 
will  outline  the  necessary  steps  which  must  be  taken  to  remove  these 
assumptions. 

The  first  assumption  is  that  of  a  linearized  elastic  wave  equation. 
This  is  widely  used  and  accepted  by  seismologists  and  is  justified  by  the 
small  amplitudes  of  differential  strains  that  are  produced  globally  by  most 
seismic  disturbances.  The  exception  to  this  is  in  the  near  vicinity  of  a  large, 
non-linear  seismic  source  such  as  an  earthquake  or  a  nuclear  explosion. 
Even  in  these  cases  we  can  consider  the  elastic  wave  field  to  be  linear 
beyond  some  volume  which  encloses  the  non-linear  source  region  and 
represent  the  source  by  an  equivalent  linear  point  source  (e.g.  Bache  and 
Harkrider  (1976)).  In  this  case  our  solution  will  be  invalid  within  the  non¬ 
linear  region  but  will  be  valid  outside  that  region. 


The  second  assumption  is  an  earth  structure  which  is  laterally 
homogeneous.  This  is  the  most  restrictive  assumption  and  considering  the 
attention  that  seismogram  synthesis  for  laterally  homogeneous  structures 
has  received  in  the  past,  some  justification  is  deserved  at  this  point.  I  justi¬ 
fied  the  laterally  homogeneous  assumption  based  upon  the  following  three 
premises: 

1.  A  starting  point  for  the  solution  of  certain  laterally  inhomogene¬ 
ous  problems  is  an  accurate,  efficient  and  complete  solution  to  the 
laterally  homogeneous  problem. 

2.  Although  much  work  has  been  directed  towards  the  laterally 
homogeneous  problem,  an  accurate  efficient  and  complete  solution 
is  yet  to  be  realized. 

3.  Spectral  solutions  of  the  elastic  wave  equation  promise  the  most 
efficient  solution  to  the  problem,  at  least  when  compared  to  the 
methods  used  presently. 

The  second  premise  can  be  justified  by  a  quick  review  of  existing  methods 
for  synthesizing  seismograms  for  laterally  homogeneous  structures.  These 
methods  fall  into  two  general  categories:  ray  theoretical  methods  and  com¬ 
plete  solutions  methods.  Ray  theoretical  methods,  although  efficient  and 
accurate,  do  not,  in  general,  provide  complete  solutions  (Hartzell  and  Helm- 
berger  (1982)).  The  complete  solution  methods  are  all  closely  related  to  the 
original  reflectivity  method  (Fuchs  and  Muller,  (1972))  and  suffer  from 
being  rather  inefficient.  This  is  due  to  the  direct  numerical  integration 
approach  used  by  all  of  these  methods  to  compute  the  inverse  Hankel 
transform.  Discrete  spectral  representations  of  the  solution,  however,  allow 
for  the  inverse  Hankel  transform  to  be  evaluated  analytically  in  terms  of  a 
residue  sum  which,  in  addition  to  avoiding  the  accuracy  problems  associated 


with  direct  numerical  integration,  also  provides  a  more  efficient  solution  as  1 
will  show  later. 

The  third  assumption  is  the  flat  earth  approximation.  This  is  not 
at  all  restrictive  and  can  be  easily  justified  for  several  reasons.  First  of  all, 
as  will  be  shown  later,  the  major  applicability  of  the  spectral  method  will  be 
for  near  field  problems  where  the  flat  earth  approximation  is  quite  good  on 
its  own.  If  the  source  receiver  distances  become  large,  a  flattening  correc¬ 
tion  can  be  applied  to  the  vertical  structural  dependence  which  will  give 
very  good  solutions  out  to  teleseismic  distances.  It  should  be  pointed  out 
that  all  ray  theoretical  solutions  end  up  malting  the  flat  earth  approxima¬ 
tion  implicitly  by  using  the  Poisson  sum  formula  to  convert  discrete  sum 
solutions  in  terms  of  spherical  wave  functions  to  continuous  integral  solu¬ 
tions  in  terms  of  cylindrical  or  Cartesian  wave  functions. 

The  fourth  assumption  is  that  the  earth  structure  will  be  made 
up  of  homogeneous  layers  connected  by  flat  plane  interfaces.  At  first  this 
may  appear  to  be  a  restrictive  assumption,  but  we  can  represent  any  arbi¬ 
trary  depth  dependence  for  some  finite  wave  length  by  specifying  a  large 
number  of  suitably  thin  layers.  The  question,  then,  is  whether  it  is  more 
efficient  and  accurate  to  use  a  large  number  of  homogeneous  layers  or  a 
smaller  number  of  inhomogeneous  layers  and  I  will  address  this  question  in 
more  detail  in  a  later  section. 

The  fifth  assumption  is  that  the  earth  structure  will  be  isotropic. 
This  restriction  can  be  reduced  to  that  of  lateral  isotropy  without  any  fun¬ 
damental  change  in  the  analytical  development  of  the  solution  although 
with  the  expense  of  significantly  increased  algebraic  complexity. 


The  sixth  assumption  is  that  the  source  of  seismic  energy  will  be 
a  spatial  point  source.  Once  again  this  is  not  a  restrictive  assumption  since 
any  distributed  source  can  be  represented  as  either  an  equivalent  point 
source  or  a  linear  superposition  of  a  large  number  of  point  sources  placed 
along  the  boundary  of  the  distributed  source. 


2.2  Differential  Equations,  Coordinate  Systems,  and  Boundary 
Conditions 

The  most  general  representation  of  the  linearized  elastic  wave 
equation  in  Cartesian  coordinates  is: 


*U.j  +  *  pu'i  (2-2.1) 

where  is  the  space  and  time  dependent  stress  tensor,  is  the  space  and 
time  dependent  applied  body  forces,  p  is  the  space  dependent  density  and  u2 
is  the  space  and  time  dependent  displacement  vector.  Integer  indices  range 
from  one  to  three  with  the  implied  summation  convention.  Partial  differen¬ 
tiation  is  indicated  by  ,  j  =  9/9xj  and  "=  92/9t2. 

The  stress  tensor  is  related  to  the  displacements  through  the  con¬ 
stitutive  tensor,  Cijkl: 


*ij  =  Cijkl  uk>| .  (2.2.2) 

Using  the  first  law  of  thermodynamics  and  assuming  adiabatic  elastic  defor¬ 
mations,  one  can  show: 


Cjjkl  _  Ck|jj  , 

and  since  the  stress  and  strain  tensors  are  symmetric, 

Cjjki  =  CjikJ  =  Cjj|k  =  Cj •, j k  . 


(2.2.3) 


(2.2.4) 


Equations  (2.2.3)  and  (2.2.4)  reduce  the  number  of  independent  components 


in  C,,n  to  21  and  this  is  the  number  of  elastic  constitutive  parameters 
which  must  be  specified  for  a  general  anisotropic  material.  Since  we  have 
assumed  a  linearized  elastic  wave  equation,  each  of  these  parameters  depend 
only  on  the  spatial  coordinates  and  not  on  the  displacement  or  any  of  its 
derivatives.  Also,  since  we  assumed  a  laterally  homogeneous  structure,  the 
elastic  constitutive  parameters  (along  with  the  density)  will  only  be  depen¬ 
dent  on  one  spatial  coordinate,  xs. 

If  we  assume  lateral  isotropy  about  the  xs  axis,  the  number  of 
independent  elastic  constitutive  parameters  reduces  from  21  to  5  and  follow¬ 
ing  the  notation  of  Takeuchi  and  Saito  (1972)  we  represent  these  parame¬ 
ters  with  the  coefficients  A,  C,  F,  L  and  N  in  expression  (2.2.2). 

=  A(uj  j  +  Uj'j)  -  2Nu2i2  +  Fu3  j  (2.2.5) 

*22  =  A(ui,i  +  u2>2)  -  2Nu, ,  +  Fuj  j 
*33  =  F(u,i  +  uw)  +  Cu33 

*23  =  L(u2i3  +  u3>2)  =  o  12 

*31  =  L(u5i|  +  Uj^)  =  <r,3 

*12  =  N(ul^  +  u2.l)  =  *21 
For  the  case  of  a  completely  isotropic  material, 

A  =  C=  A  +  2*,L=N=/i,F=A  (2.2.6) 

where  A  and  p  are  the  Lame  elastic  parameters  and  we  can  rewrite  equation 
(2.2.2)  in  the  following  familiar  form; 

*u  =  A  *,juk,k  +  n(uui  ■+  uj  i)  .  (2.2.7) 

Typically  the  elastic  parameters  are  redefined  in  terms  of  the  P  and  S  wave 


velocities  in  a  homogeneous  isotropic  material 


o  =  v/(A  +  2 n)/p  ,  0  -  s/iTfp  (2.2.8) 

For  the  case  of  a  laterally  isotropic  material  we  can  also  define  the  elastic 
parameters  in  terms  of  horizontal  and  vertical  propagating  P  and  S  wave 
velocities. 

qh  =  VAjp  ,  ov  =  , 

(2.2.9) 

0H  =  v'NTp  ,  0\  =  vt/p 

and  we  are  left  over  with  a  fifth  coefficient,  F. 

The  constitutive  relations  (2.2.5)  and  (2.2.7)  along  with  the  elas¬ 
tic  parameters  allow  us  to  model  either  a  full  elastic  material  or  an  acoustic, 
liquid  material  by  setting  0-0.  We  can  also  model  an  anelastic  material 
by  allowing  the  elastic  moduli  to  have  non-zero  imaginary  components  and 
this  will  be  addressed  in  more  detail  in  a  later  chapter. 

In  order  to  easily  represent  the  radiation  field  from  point  sources 
and  to  match  boundary  conditions  at  horizontal  layer  interfaces,  the  coordi¬ 
nates  will  be  changed  from  Cartesian  to  a  cylindrical  coordinate  system  and 
all  of  the  following  analytical  developments  will  be  done  in  cylindrical  coor¬ 
dinates.  The  cylindrical  coordinate  system  is  shown  in  Figure  2-1  along 
with  layer  numbering  conventions.  At  a  later  point  in  the  theoretical 
development  the  assumption  of  a  layered  structure  will  be  made  and  since  I 
will  be  writing  solutions  of  the  wave  equation  in  each  individual  layer  and 
then  matching  boundary  conditions  throughout  the  stack,  I  employ  both  a 
global  coordinate  system  and  a  set  of  local  coordinates,  each  relative  to  an 
individual  layer.  The  origin  of  the  global  coordinate  system  will  be  at  the 
free  surface  with  the  positive  z-axis  pointing  down.  The  origin  of  a  local 


Figure  2-1.  These  are  the  coordinate  systems  and  layer  numbering  conven 
tions  which  are  used  in  the  theoretical  development. 
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coordinate  system  will  be  at  the  top  of  the  layer,  with  the  radial  and  azimu¬ 
thal  coordinates  being  the  same  for  all  of  the  coordinate  systems. 

I  distinguish  between  the  global  and  local  vertical  coordinates  by 
using  an  unsubscripted  or  unsuperscripted  z  for  the  global  coordinate  and  a 
superscripted  for  local  coordinates  where  p  is  the  layer  index.  The 
depth  of  the  bottom  of  the  pth  layer  in  global  coordinates  is  z  =  whereas 
the  thickness  of  the  pth  layer  is  (£(p)  =  h(pl  -  h{p_1)).  This  dual 

representation  is  used  for  all  of  the  functions  of  z  as  well.  Whenever  a  func¬ 
tion  of  z  appears  without  a  layer  superscript  it  is  understood  that  the  argu¬ 
ment  will  be  in  global  coordinates  and  whenever  a  layer  superscript  does 
appear,  then  the  argument  of  the  function  will  be  in  local  coordinates. 
Thus  for  some  function,  f(z) 

f«  I  ,rt)  (2.2.10) 

By  making  the  flat  earth  and  laterally  homogeneous  assumptions,  the  elastic 
moduli  and  density  will  be  dependent  only  on  the  z  coordinate  and  the 
lateral  isotropy  assumption  will  imply  isotropy  about  the  z-axis. 

The  boundary  conditions  will  be  specified  either  in  terms  of  dis¬ 
placements  and  tractions,  or  in  terms  of  wave  field  functions  (radiation  con¬ 
dition).  The  displacements  are  denoted  by  ut(r,0.z:t),  u*(r,0,z;t),  Un(r,0,z;t), 
and  the  tractions  across  a  horizontal  plane  are  Tr(r,0.z;t),  T#(r,0,z;t),  and 
Tt(r,0,z;t).  The  boundary,  conditions  are  as  follows: 

1.  A  traction-free  surface  will  exist  at  the  top  of  the  structure  (z  =  0), 

Tr(r,*,0;t)  =  T,(rA0;t)  =  (2.2.11) 

=  Tt(r,0,O;t)  *  0  for  all  r,0,t  . 

2.  Either  a  radiation  condition  will  exist  for  a  semi-infinite  bottom  layer. 


no  sources  at  z  =  oc  , 


(2.2.12) 


or  a  plate  boundary  condition  will  exist  for  the  underside  of  the  bottom 
layer  at  z  «  H  which  can  be  specified  in  several  ways;  either 

T,(r,*,H;t)  =  T#(r  AH;t)  =  Tt(r,*,H;t)  -  0  for  all  r^,t,  (2.2.13) 

or 

ur(r,0,H;t)  =  u#(r,^,H;t)  =  ug(rAH;t)  =  0  for  all  r,0,t,  (2.2.14) 

or 

Tr(r,0,H;t)  =  T#(rAH;t)  =  ug(rAH;t)  =  0  for  all  r,*,t  .  (2.2.15) 

The  mixed  boundary  condition  given  by  (2.2.15)  insures  no  conversion 
of  P  to  S  energy  at  the  plate  bottom  or  vice  versa. 

Tractions  and  displacements  will  be  continuous  along  the  2-axis  for  all 
r,0,  and  t  as  long  as  the  elastic  moduli  are  continuous  except  at  the 
source  location. 

For  a  horizontal  interface  at  z  =  h  where  the  elastic  moduli  change 
discontinuously,  the  boundary  conditions  will  be  specified  according  to 
the  type  of  discontinuity.  For  a  solid-solid  discontinuity, 

T,M,h+;t)  =  Tr(r,*,h-;t)  , 

T#M,h*;t)  =  T,(rAh-;t), 

Tg(r,*,h4;t)  =  Tt(rAh-;t),  (2.2.16) 

ur(rAh+;t)  =  ur(rAh';t), 
u,(rAh*;t)  =  u«(r,0,h~;t), 
and  ug(r,0,h*;t)  =  ug(r,tf,h";t), 

for  all  r A  and  t,  where  h4  is  just  below  the  interface  and  h~  is  just 
above  the  interface. 
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For  a  solid-liquid  and  liquid-liquid  discontinuity 
Tr(r,*,h+;t)  -  T#(rAh+;t)  =  0 , 

Tf(r,0,h”;t)  =  T#(r,0,h~;t)  *  0, 

Tt(r,*,h+;t)  =  Tg(r,0,h";t),  (2.2.17) 

u,(r,0,h+;t)  =  ug(r,0,h‘;t), 

for  all  r,fi,  and  t. 


2.3  Separable  Solutions  of  the  Elastic  Wave  Equation  in  Cylindri¬ 
cal  Coordinates 

The  elastic  wave  equation  (2.2.1),  can  be  written  in  cylindrical 


coordinates  as  follows: 


,  .  1  3  ,  ,  1 

'“'■'f'  +  7  ar(r<,")T7-5T  + 


(2.3.1) 


I  | 

+  “  7'* 


,  .  1  d  .  .  1  d<r»e 

/>u, .  +  -  —  (r<r*)  +  7  — —  + 


.  .  1 

+  — —  +  — 
dz  r 


.  .  1  d  .  ,1 

+  7  -  <r.„)  +  -  —  +  _ 


Likewise,  the  constitutive  equations  (2.2.5)  can  also  be  written  in 


cylindrical  coordinates  for  the  case  of  lateral  isotropy. 


du,  i  du*  i  du. 

„„.A_L+(A-2N)  +  F^ 


(2.3.2) 


I 


,  0u.  1  8u«  i  0u. 

'-■(A-JN)-af  +  A  T-5T  +  7"-  +  f'5T 


P  3ur  ^  r  1  0U*  ^  1  ^  n  d“* 

<r  M  =  F  +  F  “  +  “  ur  +  C  ~ 

or  r  o8  r  8s 


1  8ut  8u, 

"•■L  7-dT  +  — 


dur  0ut 

'"‘l  "57  "a7 


0U#  1  I  0Ur 

'-’N  "a7  "  7  “*  7  "aT 


where  the  coefficients  A,  C,  L,  N,  and  F  are  functions  only  of  z.  For  the 
case  of  a  completely  isotropic  structure  equations  (2.3.2)  reduce  to  the  fol¬ 


lowing: 


.  0U.  1  1  0U*  0U. 

■  <A  +  -aT  +  A  7  +  7  ~af  +  "57 


,  i  1  du«  8ur  8u. 

,„.(A  +  2„)  7  *  7  -jj-  -37+-a7 


,  „  0U.  1  1  0U«  0U- 

-  (A  +  2„)  —  +A  _Ul+__+_ 


1  8u.  du« 

7  dT  +  -dT  =  4r‘' 


(2.3.3) 


V  tr  *  P 


8ur  8ut 
0r  +  0r 


0u#  1  1  dur 

'«'»  "57  ”  7  7  aT  * 


v  r.  ;.  -.v.v  < 
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where  once  again  A  and  p  are  functions  only  of  t. 

Equations  (2.3.2)  or  (2.3.3)  could  be  substituted  back  into  equa¬ 
tions  (2.3.1)  to  eliminate  the  stress  variables  and  produce  a  system  of  three 
coupled,  second  order,  partial  differentia]  equations  in  the  three  displace¬ 
ment  unknowns,  ur,  u#1  and  u,.  The  next  step,  normally,  would  be  to  define 
three  potential  functions  of  the  displacements  in  order  to  decouple  the  sys¬ 
tem  of  equations  (2.3.1).  The  P  wave  dilatational  potential  and  the  SV  and 
SH  wave  shear  potential  functions  will  only  decouple  equations  (2.3.1)  when 
the  structure  is  completely  homogeneous  and,  in  the  case  of  an  arbitrarily 
inhomogeneous  material,  the  representation  of  equations  (2.3.1)  in  terms  of 
the  potential  functions  would  be  completely  coupled.  Ben-Menahem  and 
Singh  (1981)  have  shown,  however,  how  separable  solutions  of  the  elastic 
wave  equation  in  cylindrical  coordinates  for  the  case  of  an  arbitrary  verti¬ 
cally  inhomogeneous  structure  can  be  represented  in  terms  of  vector  cylindr¬ 
ical  harmonics.  When  these  separable  functions  are  substituted  into  equa¬ 
tions  (2.3.1),  the  t,  r  and  9  dependence  are  eliminated  and  we  are  left  with 
a  set  of  coupled  ordinary  differential  equations  (ODE)  in  the  depth  variable, 
z. 

I  will  thus  assume  the  following  solutions  for  the  displacements, 
u  =  (u„u,,ut). 

u  =  t"'  |up(u?,k,m,z)  PiT(r,0)  -  (2.3.4) 

+  uB(w,k,m,*)  Bkm(r ,9)  +  uc(u>,k,m,t)  Ckm(r,0)J 

where  w  is  the  constant  angular  frequency, 

k  is  the  constant  horizontal  wavenumber, 
m  is  the  constant  integer  azimuthal  order, 


.*  'A-  '?*  '*  j 


and  P,  B  and  C  are  the  vector  cylindrical  harmonics  and  are  given  by 

P“M)  -  Yk“(r^) 

*  (*■8^) +,,tt  Jf)  Y"(,'*)  <23-5> 

c"(r’>>-(*’Tr^--8R|Y"(r’')' 

where  er,  %t  and  e,  are  unit  vectors  in  the  r,  9  and  z  directions  and  Yk“(r,0) 
is  the  scalar  horizontal  wavefunction  and  is 

Yk“M)  =  JB(kr)  eto'  .  (2.3.6) 

Jm(kr)  is  the  integer  order  Bessel  function  of  the  first  kind. 

In  order  to  satisfy  the  boundary  conditions  it  will  be  necessary  to 
compute  the  tractions  across  the  z— constant  plane  as  well  as  the  displace¬ 
ments.  Separable  solutions  for  the  tractions  in  terms  of  vector  cylindrical 
harmonics  will  also  be  used. 

T  =  twt  (Tp(W,k,m,*)  P„m(r,*)  +  (2.3.7) 

+  TB(u>,k,m,z)  Bkm(r ,9)  +  Tc(u>,k,m,z)  Ckm(r,0)J 

where  T  =  (Tr,T#,Tt)  and 
T  r  = 

T#  *  o  H 

T,  = 

Finally,  I  will  represent  the  body  forces,  pt  =  {pir,  p(t,  p[t),  in  the  same 
manner  as  the  displacements  and  tractions. 


(2.3.8) 


(fP(wk,m,z)  P*(r,0)  4 

4  fB(u>,k,m,*)  B*(r,8)  4  fc(w,k,m,z)  Ck“(r,8)j 

If  we  substitute  equations  (2.3.4),  (2.3.7),  and  (2.3.8)  into  equa¬ 
tions  (2.3.1)  and  (2.3.2)  and  use  relations  (2.3.5)  and  (2.3.6)  along  with 
various  properties  of  integer  order  Bessel  functions,  the  t,  r  and  6  depen¬ 
dences  can  be  factored  out  leaving  a  system  of  six  coupled  ordinary  differen¬ 
tial  equations  of  z,  in  the  dependent  variables  uP,  uB,  uc,  Tp,  TB,  and  Tc. 
The  resulting  z-dependent  ODEs  can  be  written  as  follows: 


Once  gain  the  elastic  moduli  A,  C,  L,  N,  and  F  are  arbitrary  functions  of  z 
and  at  this  point  we  have  made  no  approximations  or  assumptions  other 
than  linearity,  laterally  homogeneous  structure  and  laterally  isotropic  struc- 
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Because  of  the  separable  form  of  the  solutions  expressed  by  equa¬ 
tions  (2.3.4)  and  (2.3.7),  equations  (2.3.7)  must  hold  at  a  fixed  u>,  k  and  m 
for  all  t,  r  and  6.  Thus  the  boundary  conditions  expressed  by  equations 
(2.2.11)  through  (2.2.17)  will  be  insured  by  applying  them  to  the  depth 
dependent  spectral  functions  in  equations  (2.3.9). 

2.4  Solutions  of  the  Depth  Dependent  ODEs:  The  Propagator 
Matrix 

We  can  write  equations  (2.3.9)  in  a  more  convenient  form  as  fol¬ 
lows: 


(y(*)>  =  Ms)]  {y(*)>  -  M*)}  ,  (2-4.1) 

where  (y(z)}  is  the  six  component  displacement-stress  vector,  and  is  given 
by 


7i(*) 

Up 

yz(*) 

»B 

y*(*) 

Tp 

y«(*) 

T„ 

ys(*) 

uc 

y«(*) 

Tc 

i 

{w(z)}  is  the  six  component  forcing  function  vector  and  for  a  time  and  space 
distribution  of  simple  body  forces  is  given  by 


wi(*) 

0 

w2(z) 

0 

*s(*) 

fp 

wa(*) 

^B 

*s(*) 

0 

*s(*) 

fc 
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and  (U(x)|  is  a  six  by  six  matrix  whose  element  are  functions  only  of  w,k 
and  depth  dependent  elastic  moduli. 

It  is  obvious  from  equations  (2.3.0)  that  we  can  partition  equa¬ 


tions  (2.4.1)  as  follows, 


-r  =  Iru(01  {ny(*)>  -  {**(*)} . 

dz 


-j-  (iM*)} -  IlU(*)1  -  (lw(*)>  » 

dz 


Where 


(ay(*)}  - 


y»(*) 

yx(*) 

y,(«)  '<*wW) 

y«(*) 


Wj(l) 

wr(*) 

Wj(*)  * 

w«(«) 


IrU(x)]  = 


0  k*  - 


<Ly(*)> 


(ys(*)l  /ws(*)l 

■  (>«(•)( Aw(,)>  ‘  W«)|  ■ 


1lU(Z)]  =  (k^HV-u ,2p)  0 


(2.4.4) 


(2.4.5) 


(2.4.6) 


(2.4.7) 
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and  where  the  R  subscript  denotes  Rayleigh  or  P-SV  wave  motion,  the  L 
subscript  denotes  Love  or  SH  wave  motion,  and  the  elastic  moduli  A,  C,  L 
and  N  have  been  replaced  by  the  P  and  S  wave  vertical  and  horizontal  velo¬ 
cities,  ov,  0\,  aH,  and  fin  from  equations  (2.2.9).  We  can  easily  determine 
[rU]  and  [lU]  for  a  completely  isotropic  structure  by  setting  ov  =  aH  =  a, 
£v  =  =  &  and  F  =  A. 

Consider  the  homogeneous  or  unforced  part  of  equation  (2.4.1), 

namely, 

—  fr>  -  iU]  b)  ■  (2.4.8) 

Gilbert  and  Backus  (1966)  studied  solutions  to  this  equation  and  popular- 

« 

ized  and  generalized  the  propagator  matrix  method  which  was  first 
employed  by  Thomson  (1950)  and  Haskell  (1953)  for  the  case  of  a  plane 
homogeneous  layered  structure.  In  order  to  solve  (2.4.8)  the  matrizant  of 
[U]  is  defined  as  follows, 

K 

|A(.,.„)]  *  [I)+  f|U(tl)]d(,  T  (2.4.9) 

*9 

l  fl 

+  f|U(fi)l  JT  |U(f,)l  dt,dr,  +  •  •  • 

where  |I]  is  the  unit  matrix.  Differentiating  the  matrizant  with  respect  to  z 
gives, 

t 

±  |A(i,i0)J  -  |U{i)]  +  |U(.))  J|U(?,)i  dSl  +  •  •  • 

or 

4-  |A(.,*o)l  -  |UW)  |A(«*)!  ■  (2.4.10) 


If  we  post  multiply  (2.4.10)  by  {?(><>)}  we  arrive  at, 

|A(»a)1  M»o)>  -  |U(.)1  [A(.,»o)l  b-<«o)>  ,  (2.4.11) 

and  comparing  this  to  equation  (2.4.8)  it  follows  that, 

(y(t)}  =  (A(z,z0)]  (y(*0)}  •  (2-4.12) 

When  we  include  the  forcing  function  vector  {w},  equation  (2.4.12)  becomes 

S 

(y(*)>  =  [A(«,*o)l  (y(*o)}  -  (A(*,f )1  {w(f )}  df  .  (2.4.13) 

The  matriz&nt  of  [U]  given  by  equation  (2.4.9)  allows  us  to  pro¬ 
pagate  the  solution  for  {y}  from  some  initial  depth,  Iq,  to  some  other  depth, 
z,  and  thus  the  matrizant  is  usually  called  the  propagator  matrix.  One 
obvious  property  of  the  propagator  matrix  which  must  be  true  in  order  for 
relation  (2.4.12)  to  be  valid  is  that, 

{A(z,z)]  =  [I]  .  (2.4.14) 

Another  important  property  of  the  propagator  matrix  is  that  it  is  only 
dependent  on  w,  k,  z,  and  the  depth  dependent  elastic  moduli,  and  is 
independent  of  r,  9  and  the  azimuthal  order  m.  The  major  difficulty,  then 
in  determining  spectral  solutions  of  the  elastic  wave  equation  for  a  laterally 
homogeneous  structure,  is  in  computing  the  propagator  matrix. 

Consider  the  {y}  vector  at  three  depths,  z0,  ij,  and  z2,  such  that 
z2  >  z,  >  z0.  From  equation  (2.4.13), 

(y(*i)}  =  [ a(*j,*o)]  {y(*o)>  -  /  [A(*i,f )i  M? )) 

*• 

and 

*» 

{y(*x)>  =  1a(*2,zi)}  (y(*i)>  -  f  !A(*j,f)J  {w(f )}  df , 
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(y(*a)}  *  (A(*a,*i)]  (A(*i4o))  (y(*o)} 


-  IA(*j,*i)]  f  [A(*„f )]  {w(f )}  df 


-  f  |A(*2.f)]  Mf))  df  i 


but  since 


{y(*2)>  *  (A(zj,*0)]  (y(*o)}  -  /  (A(*2,f)]  Mf))  df  , 


iA(*2>*o)]  -  !A(*2»*l)]  |A(*i,*o)]  • 


(2.4.15) 


Equation  (2.4.15)  expresses  an  important  property  of  the  propagator  matrix 
which  is  especially  useful  when  dealing  with  layered  structures.  The  depth 
restriction,  z2  >  zt  >  z0,  is  actually  not  necessary  and  as  a  corollary  to 


(2.4.15), 


[A(z0,z0)]  =  [1]  =  [A(z0,Z|)J  (A(Zj,Zo);  , 


|A(*o.*i)j  -  [A(*i,*0)]"1  • 


(2.4.16) 


Thus  the  upward  propagator  matrix  between  two  depths  is  the  inverse  of 
the  downward  propagator  matrix  between  the  same  two  depths. 

For  the  case  of  a  general  heterogeneous  structure  with  depth, 
equation  (2.4.9)  cannot  be.  solved  analytically  and  the  only  method 
presently  available  is  that  of  numerical  integration.  In  practice  equation 
(2.4.1)  is  integrated  with  respect  to  z  so  that, 

ft  ft 

M«))  -  M«o»  +  f  |U(f )]  Well  if  -  ^  Mr ))  dr  . 


Wa;v.v;,;v:.\v.v.w.  V.  W  v.  v.'.v 


(2.4.17) 

The  six  by  six  propagator  matrix  can  be  determined  by  ignoring  the  forcing 
function  contribution  in  (2.4.17)  and  then  numerically  integrating  (U(z)| 
(y(z)}  six  times  with  each  element  of  the  starting  value  of  {y (*o)}  set  con¬ 
secutively  to  unity,  the  other  elements  being  zero.  The  resulting  6ix  values 
of  the  vector  (y(z)}  will  constitute  the  six  columns  of  (A(z,z0)].  It  should  be 
pointed  out  at  this  time  that  all  solutions  of  the  depth  dependent  ODEs  for 
an  arbitrary  heterogeneous  medium  are  in  fact  approximations  to  the  exact 
solution.  Thus  we  should  evaluate  candidate  approximations  according  to 
accuracy,  efficiency  and  ease  of  implementation. 

One  candidate  approximation  which  is  very  popular  and  rela¬ 
tively  easy  to  implement  is  to  assume  the  structure  is  made  up  of  a  number 
of  plane  layers,  each  layer  being  completely  homogeneous.  It  is  obvious 
that  any  arbitrary  depth  dependence  can  be  approximated  by  a  large 
number  of  sufficiently  thin  homogeneous  layers  welded  together.  The 
approximation  would  then  break  down  at  wavelengths  less  than  or  equal  to 
the  individual  layer  thicknesses,  and  the  approximation  could  be  made  arbi¬ 
trarily  accurate  by  decreasing  the  layer  thicknesses. 

In  order  to  see  why  the  homogeneous  layered  structure  approxi¬ 
mation  is  easy  to  implement  let  us  first  consider  the  case  of  a  completely 
homogeneous  structure.  In  this  case  [U]  is  independent  of  z  and  can  be 
taken  outside  of  the  integrals  in  (2.4.9)  which  then  gives 

|Ata,)1  -  [I]  +  (I  -  «J|U|  +  {  (I  -  *o)JiU',|V!  +  •  •  • 

-  exp|(z  -  *o)|U]j  .  (2.2.18) 

We  can  apply  Sylvester’s  formula  (Hildebrand  (1952))  to  compute  any  func- 
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tion  of  a  matrix  in  terms  of  the  eigenvalues  of  that  matrix. 

.  n  <iui  -  a.|i» 

(2"19) 

m#  k 

where  [U]  is  a  square  n  x  n  matrix  with  distinct  eigenvalues,  Ah  i  =  l,2,...,n. 
We  now  have  an  explicit  analytical  solution  for  the  propagator  matrix  in 
terms  of  the  eigenvalues  of  (z  -  z0)[U).  From  equations  (2.4.18)  and  (2.4.19) 
we  can  see  right  away  that  the  solutions  for  the  propagator  matrix  will  be 
exponential  functions  in  z  which  is  what  we  might  have  expected  for  a 
homogeneous  structure.  It  is  the  exact  solution  for  the  propagator  matrix 
given  by  (2.4.19)  along  with  the  simple  functional  dependence  with  depth 
that  makes  the  homogeneous  layered  approximation  easy  to  implement, 
relatively  efficient  and  arbitrarily  accurate.  As  for  as  I  can  determine,  there 
is  no  other  structural  approximation  for  which  an  exact  analytical  solution 
for  the  propagator  matrix  can  be  obtained. 

The  eigenvalues  of  (z  -  z0)[U]  for  SH  waves  are  easy  to  compute 
and  from  equation  (2.4.7)  are 

a  2  j  U/2 

.  --jt) 

a  2  j  '1/2 

A* --(*-«•) 

P\  P\ 

The  Love  wave  propagator  matrix  can  now  be  expressed  as  follows  for  a 
solid  laterally  isotropic  homogeneous  material. 

[LA(z,z0)]  =  exp  [(z  -  z0)  IlU;]  , 


=  exp(A,) 


(IlU)  ~ 

(A,  -  A2) 


(2.4.21) 


■  '*'  JftC  J  fi  '  i  >4  , 


■jfcVjto.afc-. 


A 


-  exp(Aj) 


(IlU|  -  a ,[!!) 


(A,  -  A2) 


IlA(*,*o)]  =  tU]  -j-L-  «inh(i(*  -  *o)  tr ) 


+  (I]  cosh  |i(z  -  »0)  i*')  , 


[LA(t,*0)J  =  [lU]  -i-  sin|(z  -  *0)  L*/) 


+  [I]  cos((*  -  *o)  tv)  , 


(2.4.22) 


where  \y  is  the  Love  wave  vertical  wavenumber  and  is  given  by. 


.(  0 H  .2 

»*/*-!  — -  k  -  — - 


(2.4.23) 


I  have  chosen  the  vertical  wavenumber  by  analogy  with  the  horizontal 


wavenumber  so  that  when  \y  is  real,  propagating  solutions  occur  in  the  z- 


direction  and  when  \y  is  imaginary,  inhomogeneous  or  evanescent  solutions 


exist.  The  propagator  matrix,  however,  remains  real  for  all  real  values  of  w. 


k,  0H,  dv  and  P-  Also,  even  though  Lv  is  a  dual  valued  function  due  to  the 
square  root,  the  propagator  matrix  itself  is  single  valued.  This  is  because 
both  values  of  Lv  are  included  in  the  propagator  matrix  which  can  be  more 
easily  seen  in  equation  (2.4.21)  where  Aj  =  +  i(z  -  z0)  \y  and 


A*  «  -  i(z  -  s0)  Li/. 

Consider  the  transformation  matrix,  |B],  which  diagonalizes  (Uj. 


A,  0  •  0 


0  A2  0 

IB]'1  tU]|Bj  =  =  [A] 


(2.4.24) 


0  0  A, 
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Obviously,  the  diagonalized  version  of  [U]  will  consist  of  diagonal  elements 
equal  to  the  eigenvalues  of  {U]  and  the  columns  of  [B]  will  be  equal  to  the  n 
eigenvectors  of  (U].  Let  us  use  [B]  to  transform  the  stress-displacement  vec¬ 
tor,  {y},  also. 

M  -  IBp1  M  , 

or 

M  -  [B]  <v>  (2.4.25) 

Substituting  into  equation  (2.4.8)  and  assuming  a  completely  homogeneous 
material, 

IB!  {*)  *  |U]|B]{v>  , 

A  {»}  -  |B]->  |U]|B]<*> , 
or 

±  {v}  -  1*1  (.}  .  (2.4.26) 

defining  the  propagator  matrix  for  the  transformed  vector,  {v},  as  [V(z,i0)j 
and  applying  Sylvester’s  formula  we  get, 

M*))  =  (V(*,*0)]  (v(s0)}  ,  (2.4.27) 

|V(*4o)j  -  "  *o)(A]j  , 


|V(*,*o)l 


exp((s  -  *0)A,) 


exp((*  -  *o)Az)  0 
0 


exp((*  -  *0)A.) 


(2.4.28) 
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This  transformed  version  of  the  stress-displacement  vector  has 
been  used  extensively  by  Kennett,  Kerry  and  Woodhouse  (1078)  and  Ken- 
nett  and  Kerry  (1979)  who  show  how  complete,  numerically  stable  solutions 
to  the  elastic  wave  equation  for  plane  homogeneous  layered  structures  can 
be  expressed  in  terms  of  a  generalized  ray  expansion.  They  refer  to  the  vec¬ 
tor,  {v},  as  the  wave  or  wavefield  vector  and  its  elements  are  decoupled 
from  each  other  in  a  homogeneous  material  due  to  the  diagonal  form  of  the 
wavefield  propagator  matrix,  (V).  The  six  elements  of  {v}  are,  in  fact,  the 
w,  k  and  z  dependent  factors  in  the  P,  SV  and  SH  potential  functions  for 
upward  and  downward  propagating  waves.  Thus  the  elements  of  the  wave- 
field  vector  are  the  w,  k  and  z  dependent  complex  amplitudes  of  P,  SV  and 
SH  rays  travelling  either  upward  or  downward  through  the  homogeneous 
material.  The  matrix  [B],  then,  transforms  the  z-dependent  solutions  from 
a  ray  representation  to  a  stress-displacement  representation.  Using  [B]  we 
can  easily  relate  the  stress-displacement  propagator  matrix  directly  to  the 
wave  field  propagator  matrix. 

[A(i,*0)i  -  [B]  |V(z,z0)]  |B;-»  (2.4.29) 

The  (B)  matrix  also  proves  useful  in  computing  reflection  and 
transmission  coefficients  across  layer  interfaces.  Because  of  the  usefulness  of 
(B],  1  will  compute  the  P-SV  propagator  matrix  using  (2.4.29)  instead  of 
Sylvester’s  formula.  The  first  step  is  to  compute  the  eigenvalues  of  [RU] 
given  by  equation  (2.4.6).  This  reduces  to  solving  for  the  roots  of  the  fol¬ 
lowing  characteristic  equation: 

bA4  +  RAJ  Ri/0J  t  ri + 

+  (r^.W  +  kWo  -  »?2))  *  0 


JL  (o’(1  -  +  l) 


(2.4.30) 


where  RA  is  an  eigenvalue  of  [gU],  and 


®  *  ®v  *0  m  fi\  «*l  “  ■'  t(  “  —  -  ®v  +  (2.4.31) 

®v  P 


k.  »  —  ,k, 
®v 


R".2 


V 

0\  ’ 

k»  -  k2  »R*//  *  kj»  -  k2 


Ri/a  and  are  P  and  S  wave  vertical  wavenumbers  and  when  they  are 
real,  the  z-dependent  solutions  are  propagating  and  when  they  are  ima¬ 
ginary,  the  z-dependent  solutions  are  evanescent  as  with  \y.  For  a  com¬ 
pletely  isotropic  material, 


9  *  Ut  -  0, 

and  we  can  write  equation  (2.4.30)  as  follows: 

rA4  +  RAJ  (ri/J  +  Ri '})  +  Ri -  0  ,  (2.4.32) 

which  can  be  factored  as, 

(rA2  +  R*'2]  [rA2  +  r^)J  *  0  , 

so  that 

rA2  -  -  r*' 2  (2.4.33) 

or 

rA2  *  -  . 

This  then  gives  the  following  four  solutions  for  RA  which  corresponds  to 
upward  and  downward  travelling  P  and  SV  waves. 

rAj  *  +  it/a  m  +  iyi/  -  kJ  “  4  rA#  (2.4.34) 

rA2  «  -  ii/m  «  -  iyij  -  ka  -  -  RA. 


rA,  «  +  w9  «  +  iyi/  -  kJ  -  +  rA^ 


">W'j 


i  ''t  J  «  l  »  1  I 


tA4  -  -  »•>  -  -  i\4/  "  k*  *  -  R*# 

We  can  certainly  obtain  an  explicit  algebraic  solution  for  RA*  even  for  the 
anisotropic  case,  however,  the  form  of  the  solution  will  generally  be  compli¬ 
cated.  For  the  isotropic  case  the  roots  RA2,  are  always  real  and  distinct  for 
real  w,  k,  o,  0  and  p  and  as  long  as  o  *  $.  For  the  anisotropic  case  the 
roots,  RA5,  may  be  either  real  or  complex  depending  on  the  sign  of  the 
discriminant  of  equation  (2.4.30).  As  long  as  17  is  close  to  one  and  f  is  close 
to  zero  (weak  anisotropy)  we  would  expect  the  roots  to  remain  real  and  dis¬ 
tinct  and  the  resulting  wave  motion  would  correspond  approximately  to  P 
and  SV  type  wave  motion.  There  will  be  values  of  17  and  however,  for 
which  the  roots,  RAJ,  will  be  complex,  but  since  the  coefficients  in  equation 
(2.4.30)  are  always  real,  complex  values  of  RA2  will  occur  as  complex  conju¬ 
gate  pairs. 

In  order  to  compute  foB]  we  need  to  compute  the  eigenvectors  of 
jRUj.  We  will  compose  |RU]  from  the  four  eigenvectors  of  [RU]  as  follows: 

t»B|  -  [(rWi  I  {**>):  |  (kM,  i  (jb),]  ,  (2.4.3b) 

where  fob}  is  a  column  vector  with  four  components  equal  to  the  ith  eigen¬ 
vector  of  [RU],  and 


||rU]  ~  rAi[I)|  fob},  =  (0)  ,i  -  1,2, 3, 4  . 


(2.4.36) 


The  normalization  of  the  eigenvectors  is  arbitrary  and  so  we  will  chose  the 
first  component  in  each  eigenvector  to  be  unity  so  that 
rB„  =  rB|2  =  rB,j  =  rBu  -  1.  Solving  for  the  remaining  components  of  the 
eigenvectors  in  (2.4.36)  we  get 


iB,:  =  1 


(2.4.37) 
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tBji 


_1_  v9  +  aV  _  ka 

Q\$\P  a9  p 


RBJi 


_1_ 

*i 


kFu;3  ,  .  k»F 

Ov^vp  1  *  «vP 


R®«i 


1  R*^*  .  ,  3  .  RAik  F 

7“  — —  +  *A>  +  — 5 — 

i;  of  Ofp 


where  S{  - 
Since  RA2 


Mi 

<*v 


0vP 


+ 1 


,  and  i  *=  1,2, 3, 4. 


-  rA,  and  rA4  -  -  rAj,  we  can  see  that, 


rBh  *  “  R®21  »  R®24  =  ~  rBm, 


(2.4.38) 


rB«  *  -  rBjj  ,  RB34  -  ~  rB33i 
R®42  *  rB«i  ,  rB«4  «  rB4S  . 

We  can  obviously  compute  |rB)_1  and  from  equation  (2.4.29)  obtain  an 
explicit  solution  for  the  stress-displacement  propagator  matrix  for  the  aniso¬ 
tropic  case. 

At  this  point  I  am  going  to  assume  a  completely  isotropic  struc¬ 
ture.  This  assumption  is  being  made  primarily  to  simplify'  the  form  of  the 
solutions  for  (B]  and  (A]  which  will  result  in  a  computationally  efficient 
algorithm,  however,  there  is  no  fundamental  reason  why  the  anisotropic 
case  cannot  be  handled  in  the  same  manner  as  the  isotropic  case.  1  am  also 
going  to  redefine  the  stress-displacement  vector,  {y},  as, 


foyf*) 


yi(*) 

ys(*) 

y»(*)A  ’ 
y*(*)A 


(2.4.39) 
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and 


Similarly,  I  redefine  the  forcing  function  vector,  {w},  as 


{R*(y)> 


wi(*) 

w2(i) 

ws(z)A 

w4(«)/k 


,{lw(z)>  - 


(2.4.46) 


1  will  drop  the  overbar  in  the  following  development  in  order  to  simplify  the 
notation.  From  the  definition  of  the  matrix,  [rB],  and  assuming  an  isotropic 
homogeneous  structure  we  arrive  at 


IrB]  = 


i*.  -»*.  1  1 

1  1  \4>f  -\*0 

pc*(Tr- 1)  *c5(Tr-l)  pc7r/i4>0  ~pc7Y\ 4>0 

pc7l\<t>a  P*\l- 1)  pcj(7-l) 


where  we  have  changed  the  normalization  from  that  defined  by  (2.4.37) 
and, 

c  =  is  the  horizontal  phase  velocity,  (2.4.42) 

4>a  =  n/(c/q)2  -  1 , 4>fi  -  >/(c//9)s  -  1  are  the  cotangents  of  the  incidence 
angles  for  plane  propagating  P  and  S  waves, 
and  7  «  2(0/c)7. 

Note  that  va  *  k^a  and  «=  kj>8  and  that  by  redefining  the  stress- 
displacement  vector  with  equation  (2.4.39)  we.  have  eliminated  us  depen¬ 
dence  from  (rB|.  We  can  easily  compute  (rB^-1  from  (2.4.41) 


7 


I 

1 

« 


1*B]-1  -  t 


-h-'y 


(7-1) 

»*• 


PC‘ 

PC‘l9m 

1 

_  1 

7 

'pc2 

pc2i*. 

(7-*) 

1 

1 

i  9ff 

pc*\9$ 

PC2 

1 

1 

i  9f 

pc7i9f 

~  pc2 

(2.4.43) 


Using  equations  (2.4.28),  (2.4.34),  (2.4.41),  (2.4.43),  and  (2.4.29)  we  can 
solve  for  the  Rayleigh  wave  stress-displacement  propagator  matrix  for  a 
solid  material. 


rA|i(s,*o)  =  -  (7-1)  «»(*.)  +  7c°s(*,)  (2.4.44) 

rA,j(*4o)  *  -  7*ft»in(ff.)  - 

9fi 

rAij(Mo)  =  -^7  *»«(*.)  +  -57-  «“(**) 

pc  p9fi 

rAh(*,*o)  =  cos(*.)  "  “T  «»(**) 

pc*  pc 

rA21(i,*o)  =  -  sin (9a)  -  ■7#jsin(0f) 

▼*0 

RAjj(*4o)  =  7  cos(*.)  -  (7-1)  cos(tf) 


rA2s(*»*o)  —  “  rAi4(*,*o) 

rA24(*4o)  =  - J—  *»"(*.)  -r  -^7  »in(**) 

peVo  Pc 


rA„(*,Xo)  =  -  P*  »*n(^0)  -  *cV**sin(f,) 

9a 


*Am(Mo)  -  ^€*7(7-  1)  cos(0#)  -  *c*7(7-1)cos(#,) 


»Am(i4o)  *  rA»(**1o) 

11^44(140)  *  ~  rAji(*,*o) 
rA4i(«,*o)  =  ~  rAj2(*i*o) 

RA42(*»*o)  =  -  pcJ72*„«n(0  -  -pC  sin(*,) 

rA4s(*,*o)  *  -  rAi2(*,*o) 


rA44(*»*o)  =  rA22(*>^0)  » 


where 


*«  =  (l  ~  *0)  =  (*  -  *0)  k*.  ,  (2.4.45) 


and 


*/»=(*-  *0)  *>  =  (*  ~  *0)  ktf ,  , 

Once  again,  as  with  the  SH  case,  we  can  see  that  although  4>a,  va,  6p  and 
are  dual  valued,  the  elements  of  the  propagator  matrix  are  all  single  valued. 

In  order  to  compute  the  propagator  matrix  for  a  liquid  or  acous¬ 
tic  layer  let  us  renormalize  the  last  two  columns  in  [rB]  with  l/(i^)  and 
then  let  $  —  0.  We  can  see  from  (2.4.41)  that  the  last  two  columns  and  the 
last  row  of  |RB]  will  go  to  zero  and  the  4x4  [RB]  matrix  will  be  singular 
and  non-invertible.  For  an  acoustic  layer  only  the  upward  and  downward 
P-wave  solutions  will  exist  and  so  we  need  to  partition  out  a  2  x  2  non¬ 
singular  matrix  from  the  4x4  [rB]  matrix  in  order  to  determine  a  solution 
for  the  acoustic  propagator  matrix.  We  define  the  acoustic  [B]  matrix  as 
follows. 


|y#(*)J  "  ^  |va(*)j 
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(2.4.46) 


where 


UB] 


i*.  -»*. 


-pc7  -pc7 


(2.4.47) 


We  should  note  that  the  second  and  third  rows  of  [RB]  are  linearly  depen¬ 
dent  for  0  *  0  and  so  for  an  acoustic  layer, 

ys(s)  -  -pc7  y2(s)  .  (2.4.48) 

We  can  invert  ^13]  to  get 


UBI"  -  j 


-TJ-  -I/PC’ 


-l/pc1 

»*• 


(2.4.49) 


and  obtain  the  2x2  acoustic  propagator  matrix. 

(yt(*) 


UA(*,*o)]  = 


!  !aA(i,*o)1 

jyi(*o) 

(ys(*o) 

(2.4.50) 

cos(40) 

/>c2sin(0a) 

-^0sin(f0) 

-pc2cos(ff0) 

(2.4.51) 

We  can  also  rewrite  (2.4.51)  in  terms  of  a  4  x  4  matrix  using  equation 
(2.4.48). 


UA(«,«o)] 


cos(fia)  0  0 

sin(tf.) 


cos(0a)  0  0 


-  ■£—  sin (Ba)  -pc7cm($m)  0  0 


0  0 


(2.4.52) 
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We  have  already  dealt  with  the  Love  wave  propagator  matrix, 
but  we  will  write  down  the  solutions  for  [lB],  [|B]_i,  and  liAfs^o)]  for  the 
case  of  an  isotropic  homogeneous  material. 


1  1 


(2.4.53) 


IlB]-1= _ ,  (2.4.54) 

P07itfi 


IlA(*,*o)1  = 

Equations  (2.4.44),  (2.4.55),  and  (2.4.52)  give  the  stress- 

displacement  propagator  matrices  for  P-SV,  SH,  and  acoustic  materials 
which  are  completely  homogeneous  and  isotropic.  In  order  to  compute  pro¬ 
pagator  matrices  for  a  homogeneous  layered  structure  we  need  to  apply  the 

layer  interface  boundary  conditions  given  by  equations  (2.2.16)  and  (2.2.17). 

(p) 

First  we  start  by  defining  the  layer  propagator  matrix,  j  CL  ],  which  relates 
the  stress-displacement  vector  at  some  point  within  the  p‘b  homogeneous 
layer  to  the  stress-displacement  vector  at  the  top  of  the  layer. 

(p) 

frlp)(f(p))}  .  [O  (t“)|  {y“(0)l,  (2.4.56) 


cos('')  7*7 

-p074>fi%\ti{0fi)  cos [9f) 


(2.4.55) 


|a<P)(?(p))]  -  !A(.,h^-»))  ,  (2.4.57) 

where  =  x  -  h<p-1\  0  ^  f^p'  <  £*p\ 

and  f*p\  ^  and  are  shown  in  figure  1. 

When  j*p*  =  £,p)  then  the  layer  propagator  matrix  will  relate  the  stress- 


displacement  vector  at  the  bottom  of  a  layer  to  the  vector  at  the  top  of  the 


layer. 


For  a  solid  material  we  can  now  start  at  any  layer  interface,  p, 


and  propagate  the  stress-displacement  vector  to  any  other  interface,  q, 
(where  p  <  q)  by  applying  the  welded  interface  boundary  conditions 
expressed  by  equations  (2.2.16)  and  by  repeating  equation  (2.4.S6)  and  in 
doing  so  we  define  the  interlayer  propagator  matrix,  |A(q  p)), 


(y(h^)}  =  [A<q*>]  {y(b(p))}  . 


(2.4.58) 


The  interlayer  propagator  matrix  is, 


q  f  (q+P  +  1  -/) 

IA(,,P)]  =  ft  [A  (£<q+p*i-0) 


(2.4.59) 


i  «p+i 


For  a  solid-liquid,  liquid-solid,  or  liquid-liquid  interface,  the  situation  is  a 
little  more  complicated  and  I  will  cover  these  cases  in  the  next  section. 

Equations  (2.4.44),  (2.4.52),  (2.4.55),  (2.4.57),  and  (2.4.59)  give 
the  ezact  solution  for  the  stress-displacement  propagator  matrix  in  a  verti¬ 
cally  heterogeneous  material  made  up  of  plane  homogeneous  layers.  Thus 
the  approximation  made  here  is  in  the  representation  of  the  structural 
model  and  not.  in  the  solution  itself.  As  I  stated  previously,  candidate 
approximations  for  the  propagator  matrix  in  an  arbitrarily  heterogeneous 
material  with  depth  must  be  compared  on  the  basis  of  accuracy,  efficiency 
and  ease  of  implementation.  The  plane  homogeneous  layered  approxima¬ 
tion  is  exact  for  the  structural  model  it  represents  and  can  be  made  arbi¬ 
trarily  accurate  to  represent  any  structural  model  in  a  straightforward  and 
physically  interpretable  manner.  Also  this  approximation  is  relatively  effi¬ 
cient  and  easy  to  implement  due  to  the  simple  algebraic-trigonometric  form 
of  the  solution. 


“V 


1  was  able  to  compare  the  homogeneous  layered  approximation 
directly  with  the  full  wave  theory  or  Langer  approximation  which  has  been 
used  extensively  by  Paul  Richards  and  was  modified  by  Vernon  Cormier 
(1980)  to  compute  the  stress-displacement  propagator  matrix  in  a  layered 
structure  in  which  the  layer  elastic  moduli  vary  linearly  with  depth.  The 
Langer  approximation  ran  about  ten  to  twenty  times  slower  than  the  homo¬ 
geneous  layered  approximation  on  a  per  layer  basis.  Also  the  Langer 
approximation  is  an  approximation  of  the  solution  and  not  in  the  represen¬ 
tation  of  the  structure  and  it  breaks  down  when  velocity  gradients  within  a 
layer  become  large.  When  this  happens  the  structure  must  be  broken  up 
into  thin  layers  as  with  the  homogeneous  layered  approximation. 

2.5  Integral  and  Spectral  Representations  for  the  Solution  of  the 
Elastic  Wave  Equation 

Equation  (2.3.4)  represents  a  solution  of  the  elastic  wave  equation 
for  all  t,  r,  0 ,  and  z  in  terms  of  the  constant  parameters  w,  k  and  m.  The 
final  solution  will  be  some  appropriate  linear  combination  of  solutions  of 
form  (2.3.4)  spanning  the  range  of  the  parameters  v,  k,  and  m,  and  this 
appropriate  combination  will  be  determined  by  the  source  vector, 
{w(w,k,m,i)},  in  equation  (2.4.13).  In  order  to  determine  the  source  vector 
we  first  define  the  following  transforms: 

40C 

^  (f(t))  =  J  f(t)  e-'1"1  dt,  -  oo  $  w  ^  +oo  (2-3.1) 

-oe 

4  00 

f'lM)  -  —  J  f(")  «+“‘  <l".  (2-5.2) 

“00 

where  F(‘)  is  the  integral  Fourier  transform  and  /,~1(*)  is  the  inverse 
integral  Fou  r  transform, 


G  (f(#))  -  J  t{*)  «‘im#  <W,  m  -  -oo,  •  -1,0,1,  •  *  oc  (2.5.o, 
c->(f(m)).^.  £  («m)  «+”*)  ,  (2.S.) 

•"  n  •  —  00 

where  G  (•)  is  the  discrete  Fourier  transform  G_1(  )  is  the  inverse  discrete 
Fourier  transform, 

00 

B( f(r))  -  f  f(r)  Jm(kr)rdr  ,  (2.S.S) 

0  ^  k  $  -foe,  m  *  -oo,  •  •  •  ,-1,0,1,  •  •  •  ,-foo 

00 

H  -1(f(k,m))  =  |  f(k,m)  Jm(kr)kdk  ,  (2.5.6) 

H  (•)  is  the  Hankel  transform  and  H  ~1(-)  is  the  inverse  Hankel  transform. 
Note  that  these  transforms  are  normalized  so  that, 

c-‘(c(W)) 

if -1  |/f  (f(r))j  *  f(r)  ,  for  all  m. 

Ben- Men  ahem  and  Singh  (1981)  show  how  any  vector,  *(r,f),  can 
be  expanded  in  terms  of  vector  cylindrical  harmonics 

00 

xM>  -  £  -J7  ,[  Mk  [xp(k,m)  Pkm(r,»)  + 

*B(k.m)  Bkm(r.f)  +  icCk"(r,»)]  ,  (2.S.7) 


\ 


where  using  the  vector  cylindrical  harmonic  orthogonality  relations, 
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tc  ft 

XP(k,ni)  -  f  rdr  f  it  *(r,«)  ■  Pk“'(r Jt) ,  (2.5.8) 


XB(kjn)  m  jri'  f  <•»  x(r.«)  •  Bkm'(r,*) 


rir  j it  ■  C”'(r.«) . 

and  *  denotes  complex  conjugate.  We  can  express  the  body  force  vector, 
lf(te>,r,#,s),  as  (2.5.7)  and  thus  the  displacement  vector  u(u/,r,6,i),  as 


*c(k,m)  =  J 


u(w,r,tf,*)  -  £  J  kdk[up(u/,k,m,z)  PkmM)  +  (2.5.9) 


+  uB(w,k,m,i)  B^r,*)  +  uc(v,k,m,z)  Ckm(r,0)|  . 

In  order  to  obtain  the  time  domain  solution,  we  apply  the  inverse  integral 
Fourier  transform  to  u(w). 


400 

u(t,r,0,z)  =  — f  d(t/  e**"1  u(u>,r,0,z)  (2.5.10) 

2>r 

The  actual  values  of  the  stress-displacement  vector,  {y}  and  thus 
the  values  of  Up,  ufi  and  uc  will  be  functions  of  fp.  fg  and  fc  in  equation 
(2.3.8)  and  using  (2.5.8)  these  are, 


00  2»  6 

fp(u>,k,m,z)  =  J  rdr  J  dfi  pf(u>,r,0,x)  •  Pkm  (r,0)  ,  (2.5.11) 


oo  2r 


fB(w,k,m,z)  =  J  rdr  j  d 8  •  Bkm  (r,f)  , 


oo  2t 


fc(w,k,m,z)  -  J*  rdr  J*  M  pf[u,r,9,x)  •  Cf*  (r,0) 


where 


pf(u>,r,#,z)  -  J  d«  e”""1  . 

Equations  (2.5.11)  allow  us  to  express  any  space  and  time  distribution  of 
body  forces  in  terms  of  frequency  and  depth  dependent  vector  cylindrical 
harmonics.  Using  equations  (2.4.3)  and  (2.4.13)  along  with  the  expressions 
for  the  propagator  matrix  and  the  boundary  conditions  at  the  top  and  the 
bottom  of  the  structure,  we  can  compute  the  stress-displacement  vector  at 
any  depth.  Finally,  with  equations  (2.4.2),  (2.5.9)  and  (2.5.10)  we  can  com¬ 
pute  the  displacement  vector,  n,  for  all  space  and  time. 

At  this  point  I  will  make  the  following  assumption  regarding  the 
source  vector,  {w}. 

{w(w,k,m,i)}  *  6{t  -zg)  {E(u,k,m)}  (2.5.12) 

I  am  thus  restricting  the  source  to  a  horizontal  plane  at  depth  za.  We  can 
see  from  equation  (2.4.13)  that  if  zg  <  z  <  zQ,  or  if  *,  <  *q  <  *>  or  if 
zg  >  *  >  zQ,  or  if  *g  >  z0  >  z  (i.e.  zg,  is  not  between  z  and  zQ),  then  the 
integral  in  (2.4.13)  will  be  zero.  On  the  other  hand,  if  z  <  zg  <  zQ,  or 
*0  <  zg  <  z,  then 

(y (*)>  =  (A(z,z0)]  (y(*0)>  -  fA(z,zg)]  {£}  ,  (2.5.13) 

z  <  zg  <  Zq  or  *Q  <  z,<  z  . 

The  vector  function  {£}  is  called  the  source  jump  vector  since  it  causes  a 
discontinuous  jump  in  the  stress-displacement  vector.  Notice  that  this 
"jump"  condition  only  exists  when  the  source  is  restricted  to  a  horizontal 
plane  (or  a  point)  and  that  for  a  spatially  distributed  source  in  depth,  the 


stress-displacement  vector  will  remain  continuous  everywhere. 

The  problem  of  computing  the  stress-displacement  vector  for  a 
given  source  at  depth  t,  has  now  been  reduced  to  a  linear  algebra  problem. 
For  a  completely  solid  material  we  can  relate  the  stress-displacement  vector 
at  the  top  interface  to  the  stress-displacement  vector  at  the  bottom  inter¬ 
face  as  follows: 

<y(0)>  *  [A(0,*t)]  (y(*  “)} ,  (2.5.14) 

(y(H)}  «  |A(H,«.)]  (y(*,+)>  , 

)>  -  M*D>  -  <£> . 

where  i"  is  immediately  above  the  source, 

and  t,4  is  immediately  below  the  source, 

and  0  *  *T  is  the  depth  of  the  top  boundary, 

and  H  -  is  the  depth  of  the  bottom  boundary. 

We  now  need  to  apply  boundary  conditions  at  the  top  and  bottom  of  the 
structure.  These  boundary  conditions  can  be  expressed  in  terms  of  zeroing 
out  some  linear  combinations  of  the  stress-displacement  vector  and  so  we 
define  the  [E]  matrices  as, 

|tE]  {,(0))  -  {0} ,  (2.5.15) 

[bE]  (y(H)}  *  (0) , 

where  the  superscript  T  denotes  the  top  interface  and  B  denotes  the  bottom 
interface. 


The  two  {£]  matrices  will,  in  general,  not  be  square  and  may  have 
different  dimensions.  For  a  solid  structure  both  [TEj  and  [®E]  will  have  six 
columns  and  the  sum  of  the  number  of  rows  for  both  matrices  will  be  six. 
The  number  of  rows  for  each  matrix  will  be  equal  to  the  number  of  boun¬ 
dary  conditions  at  that  interface  and  in  practice  there  will  be  three  boun¬ 
dary  conditions  at  each  of  the  top  and  bottom  interfaces.  The  [E]  matrices 
are  given  below  for  the  various  boundary  conditions  expressed  by  equations 
(2.2.11)  through  (2.2.15). 

1.  A  traction  free  surface, 

0  0  1  0  0  0 

[E]  *  0  0  0  1  0  0  (2.5.16) 

[o  0  0  0  0  lj 

2.  A  rigid  surface, 

1  0  0  0  0  0 

(EJ  *  0  1  0  0  0  0  (2.5.17) 

0  0  0  0  1  0 

3.  Zero  shear  tractions  and  zero  vertical  displacement, 

1  0  0  0  0  0 

0  0  0  1  0  0  (2.5.18) 

0  0  0  0  0  1 

4.  No  upward  propagating  P  or  S  wave  radiation  (Sommerfield  radiation 
condition  for  a  bottom  half  space).  In  order  to  do  this  we  must  first 
transform  the  stress-displacement  vector  to  the  wavefield  vector  using 
the  transformation  matrix,  (B]_1.  We  then  pick  out  the  first,  third  and 
fifth  rows  of  [B]-1  to  give  for  a  solid  half  space, 


JrJl 


7 - r  — : -  w  v 

pc*  pc*i+Q 

-  _J - L  o  0  ,  (2.5.10) 

l*fi  P*7 


where  the  elastic  moduli  are  those  of  the  half  space. 


No  downward  propagating  P  or  S  wave  radiation  (Sommerfeld  radia¬ 


tion  condition  for  a  top  half  space).  In  this  case  we  pick  out  the 


second,  fourth,  and  sixth  rows  of  [B]' 


fa -ril  ~  -JL  - 

•  j  '  n 


pc7\tt 


(E 1  -  ±  7  h=H  -±-  -  J- 

2  '*0  pc2i^^  pc2 


.  (2.5.20) 


Pp7'<t> , 


As  with  the  propagator  matrix,  it  is  obvious  that  we  can  partition  the 


[E]  matrix  into  a  2  x  4  [RE|  matrix  and  a  1  x  2  J^E]  matrix. 


From  equations  (2.5.14)  and  (2.5.15)  we  arrive  at, 


rE]  (A(0,if)j  (y (*,“)}  =  (0) 


(2.5.21) 


|bE]  |A(H,.,)|  fr(.,+))  -  (0) , 


[BE]  |A(H,l,)]  {y(.-)}  -  |BEJ  lAIH,.,)]  (E> 


(2.5.22) 


we  now  define  the  (D]  matrices  as  follows, 


[tD(i)J  -  |TE]  ;a(0,i);  , 


(2.5.23) 


|bD(.)]  -  |*E)  |A(H,i)]  , 


(2.5.24) 


'W. 


(2.5.25) 


(TD(*)j  {*(*)}  *  W  **  <  *,  » 

and 

(bD(«)|  {y(*)>  *  (0}  4  >  (2.5.26) 

Equations  (2.5.25)  and  (2.5.26)  are  important  relations  and  will  constitute 
the  basis  for  numerically  stable  computations  of  the  stress-displacement  vec¬ 
tor.  We  can  now  specify  the  stress-displacement  vector  immediately  above 
the  source  as, 

|TD(i,)]  {>(.-))  -  (0)  ,  (2.5.27) 

[BD(01  <y  («,“))  -  lBD(.,M  <z>  . 

This  gives  us  six  equations  in  the  six  unknowns,  {y(t~ )},  so  we  can  solve  for 
(y(if-)}.  We  can  then  compute  the  stress-displacement  vector  at  any  other 
depth  by  using  the  propagator  matrix 


{0} 

PM  M*. 

(2.5.28) 

iB(«,)l  {£) . 

<y(*)}  *  |A(s,*,))  {y(s,“)}  ,*<*,, 

(2.5.29) 

<,(.)}  -  |a(m,)]  (fr(«.-))  -  <£>)  .*  >  ». . 

(2.5.30) 

where  |D(*f)]  is  a  six  by  six  matrix  composed  of  |TD(*#)]  and  [BD(*g)]  as,  . 
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Once  Again  |D(zf)J  can  be  partitioned  into  a  four  by  four  Rayleigh  matrix 
and  a  two  by  two  Love  matrix. 

So  far,  in  order  to  compute  the  various  [D]  matrices,  we  have 
assumed  solid-solid  welded  interfaces  and  this  approach  must  be  modified 
somewhat  to  handle  acoustic  layers.  Of  course.  Love  waves  will  be  com¬ 
pletely  blocked  at  a  solid-liquid  interface  and  so  we  will  only  need  to 
address  the  P-SV  problem.  Let  us  first  consider  the  case  of  the  P-SV 
stress-displacement  vector  being  propagated  upward  through  a  solid  to 
liquid  interface.  In  this  case  we  will  denote  the  stress-displacement  vector 
in  the  solid  material  immediately  below  the  interface  as  {Ry}  and  the 
stress-displacement  vector  in  the  liquid  material  immediately  above  the 
interface  as  {Ay}.  From  equations  (2.2.17),  the  boundary  conditions  at  the 
interface  require  the  following: 

Ryl  =  A*1  ’  (2.5.32) 

R*S  =  A*3  ’ 

R>4  =  0  ’ 

A*4  =  0  ’ 

In  this  case  the  shear  displacement  will  generally  be  discontinuous.  Within 
the  solid  layer  there  are  generally  four  linearly  independent  components  of 
the  stress-displacement  vector  with  a  four  column  (and  usually  two  rows) 
[R  D]  matrix.  Within  the  liquid  layer  there  are  two  linearly  independent 
components  of  the  stress-displacement  vector  with  a  two  column  (and  usu¬ 
ally  one  row)  |  ®D]  matrix.  The  problem  then  is  to  apply  the  boundary  con¬ 
ditions  given  by  equations  (2.5.32)  to  determine  the  elements  of  |®Dj  from 
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the  elements  of  [^D]  at  the  interface.  This  is  a  straightforward  problem 
and  the  results  are  as  follows: 

PD1  *  R  D11  KD2i  ~  R DI2  R  °21  *  (2.5.33) 

A  °2  *  ®D1S  rBD„  -  rBD12  r  Dm  , 

where 


0 


Ia»1 


B 

A 


®*l 


*1 


(2.5.34) 


For  the  case  of  the  stress-displacement  vector  being  propagated  upward 
through  a  liquid  to  solid  interface,  the  same  boundary  conditions  apply  and 
the  resulting  [R  D]  matrix  is  as  follows: 


D 


l  * 


IS 


2  ’ 


(2.5.35) 


24 


1 


D 


14 


Bt 


Bi 


R  21  R  22 
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The  results  for  downward  propagating  [TD]  matrices  are  identical. 

Given  the  stress-displacement  vector,  we  can  write  the  integral 
equations  for  the  displacement  vector  given  by  (2.5.9)  as  Rayleigh  and  Love 
wave  components. 


(2.5.36) 


Rur(w,r,#,«)  «  ~  E  f  Jm(kr)  "  Jm+l(kr) 

»■  -  00  D 


e“B#  kdk 


Ru#(u>,r,*,s)  -  J  y2(u>,k,m,*;s#)  Jm(kr)  ^  tm*  kdk 

Ruf(w,r^f*)  -  ~  J  y1(w,k,m,*;*i)  Jm(kr)  eMn#  kdk 


+  00  °°  /  • 
Lur(w,r^,*)  -  J  y5(u;,k,m,t;*#)  Jm(kr)  1-^2- 


(2.5.37) 
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Lu,(u;,r,0,t)  *  -j7  S  f  y5(w,k,m,*;*,)  -  -j£-  Jm(kr)  +  Jm+1(kr)  e*"1'  kdk 


Lu|(wfr1#,s)  *  0  , 

where  >'|,  y3  and  y5  come  from  equations  (2.5.28),  (2.5.29)  and  (2.5.30). 

Equations  (2.5.36)  and  (2.5.37)  are  the  basis  for  the  various 
numerical  integration  approaches  which  I  refer  to  collectively  as  the  reflec¬ 
tivity  method.  First  popularized  by  Fuchs  and  Muller  (1971)  this  direct 
integration  method  has  been  modified  and  expanded  by  Kind  (1978),  Ken- 
nett  and  Kerry  (1979),  who  eliminated  certain  numerical  instabilities,  Cor¬ 
mier  (1980),  who  applied  the  Langer  approximation  to  model  inhomogene¬ 
ous  layers  and  deformed  the  contour  of  integration  to  avoid  singularities  in 
the  integrand  function  and  Bouchon  (1981),  who  established  a  spatial  sam¬ 
pling  theorem  with  respect  to  the  Hankel  transforms  and  applied  this 
theorem  in  a  discrete  wavenumber  method  for  evaluating  the  wavenumber 
integrals.  All  of  the  reflectivity  methods  have  in  common  the  direct  numer¬ 
ical  integration  of  the  wavenumber  (or  slowness)  integrals  in  equations 
(2.5.36)  and  (2.5.37). 
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Another  approach  for  evaluating  the  wavenumber  integrals  is  to 
deform  the  contour  of  integration  in  the  complex  wavenumber  plane  so  as 
to  encircle  the  singularities  of  the  integrand  function  and  then  apply  the 
residue  theorem.  The  integrals  given  by  equations  (2.5.36)  and  (2.5.37)  are 
not  amenable  to  this  since  the  Bessel  functions  blow  up  as  |  k  |  —  oo,  how¬ 
ever,  Lapwood  (1949)  and  more  recently  Hudson  (1969)  have  shown  how 
the  Bessel  functions  can  be  changed  to  Hankel  functions  of  the  second  kind 
by  extending  the  contour  of  integration  to  -oo. 

OO  *00 

£  f(k,m)  J,„(kr)  kdk  -  1  /  f(k,m)  HW(kr)  kdk  ,  (2.5.38) 

where  f(k,n)  *=  (-l)n+1  f(~k,n), 

and  H  (2)  is  the  integer  order  Hankel  function  of  the  second  kind. 

The  Hankel  functions  go  to  zero  as  |  k  |  -*  oo  and  lm(k)  <  0,  so 
the  contour  of  integration  can  be  closed  by  including  a  semicircular  arc  at 
infinity  in  the  lower  half  of  the  complex  wavenumber  plane.  We  now  need 
to  consider  the  locations  and  characteristics  of  the  singularities  of  the 
integrand  functions. 

The  singularities  of  the  Hankel  functions  are  well  known  and  so 
we  turn  our  attention  to  the  singularities  of  the  stress-displacement  vector, 
(y),  as  a  function  of  complex  wavenumber.  First  of  all  we  will  address  the 
question  of  when  {y}  is  a  multivalued  function  of  wavenumber  with  atten¬ 
dant  branch  points  and  branch  cuts.  The  only  multivalued  functions  to 
appear  in  the  propagator  matrices  or  boundary  condition  ([E])  matrices  are 
the  dual  valued  vertical  wavenumber  functions,  i/Q  and  (as  and  ^). 
In  the  case  of  the  propagator  matrices,  these  functions  always  appear  either 
as  arguments  of  even  functions  (e.g.  cos((z  -  z0)i/a))  or  in  products  or  quo- 
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tients  (e.g.  va  sin((*  -  *0)  vfl),  tin((t  -  *0)  */0)A'0)  «uch  that  the  result  is 
single  valued,  and  thus  the  propagator  matrices  are  single  valued.  The 
boundary  condition  matrices,  [^E]  and  jBE],  however,  may  be  either  single 
valued  or  multivalued  depending  on  the  type  of  boundary  condition. 

For  conditions  where  ail  incident  seismic  energy  is  reflected  for  all 
wavenumbers  (equations  (2.5.16),  (2.5.17),  and  (2.5.18)),  the  [E]  matrix  is 
single  valued  and  if  both  [TEj  and  [BEj  are  determined  by  one  of  these  con¬ 
ditions  then  we  will  refer  to  this  as  the  plate  problem.  For  the  plate  prob¬ 
lem  the  various  [Dj  matrices  will  also  be  single  valued,  as  can  be  seen  from 
equations  (2.5.23),  (2.5.24)  and  (2.5.31)  and  since  the  source  jump  vector, 
{E},  is  always  single  valued,  the  stress-displacement  vector  will  also  be  sin¬ 
gle  valued  as  can  be  seen  from  equations  (2.5.28)  to  (2.5.30).  Thus  the  con¬ 
tour  of  integration  can  encircle  the  lower  half  of  the  complex  wavenumber 
plane  without  being  required  to  detour  around  any  branch  cuts  or  branch 
points. 

The  case  of  most  interest  in  seismology  is  what  we  will  refer  to  as 
the  half  space  problem,  that  is,  a  reflectivity  boundary  condition  at  the  top 
and  a  radiation  boundary  condition  at  the  bottom  of  the  structure.  The 
free  surface  boundary  condition  will  be  applied  at  the  top  of  the  structure 
and  thus  [TE]  will  be  single  valued  and  given  by  equation  (2.5.16).  The 
radiation  condition  given  by  equation  (2.5.19)  will  specify  the  boundary 
condition  matrix  at  the  bottom  JBE],  but  in  this  case  the  matrix  will  be  four 
valued  for  P-SV  waves  due  to  the  i/a  and  functions  (in  the  form  of  4>q 
and  ^),  and  two  valued  for  SH  waves  due  to  (^).  We  can  see  then 
that  the  P-SV  stress-displacement  vector  will  be  four  valued  with  two 
branch  points  at  i/Q  =  0  end  =  0  and  two  branch  cuts  emanating  from 
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these  branch  points  and  the  SH  stress-displacement  vector  will  be  two 
valued  with  one  branch  point  at  *  0  and  one  associated  branch  cut.  The 
contour  of  integration  must  be  deformed  around  these  branch  cuts  and 
branch  points  in  order  to  stay  on  an  analytic  path  and  the  resulting  branch 
cut  integral  contributions  are  well  known  and  physically  attributable  to  the 
energy  which  "leaks"  away  into  the  bottom  half  space  (Gilbert  (1964)).  We 
have  been  somewhat  remiss  in  this  analysis  since  va  and  v ^  are  functions  of 

k2  and  so  for  every  branch  point  at  +k  there  is  another  at  -k.  Conse¬ 
quently,  there  are  actually  four  branch  points  and  branch  cuts  for  the  P-SV 
case  and  two  branch  points  and  branch  cuts  for  the  SH  case,  however  since 
the  integration  contour  circles  only  half  of  the  complex  wavenumber  plane, 
only  two  branch  cut  integrals  will  occur  for  the  P-SV  case  and  one  for  the 
SH  case.  As  we  will  see  this  symmetry  will  also  be  characteristic  of  the 
poles  of  {y}  as  well  as  the  branch  points  and  branch  cuts. 

The  remaining  singularities  of  the  stress-displacement  vector  are 
the  Rayleigh  and  Love  poles  which  occur  at  values  of  u  and  k  for  which  the 
[D]  matrix  in  equation  (2.5.28)  is  singular.  For  a  fixed  frequency  these 
poles  will  occur  at  discrete  wavenumbers,  however  in  the  (u\k)  space  these 
poles  form  continuous  functions  of  w  and  k  which  are  commonly  called 
dispersion  curves.  Thus,  in  order  to  locate  the  poles,  we  can  either  fix  fre¬ 
quency  and  look  for  discrete  poles  as  a  function  of  wavenumber  (or  phase 
velocity,  slowness,  etc.),  or  fix  wavenumber  and  look  for  discrete  poles  as  a 
function  of  frequency.  Whenever  the  (D)  matrix  is  singular,  we  can  write 
equation  (2.5.28)  as  follows: 

[D(Jf  (n,w),i#)]  (£(n,u>,ig)}  =  {0}  ,  (2.5.39) 


where  K  (n,u>)  is  the  nth  value  of  wavenumber  which  at  frequency  u  for 
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(2.5.39)  is  true  given  {E  (i§)}  ^  {0}.  Note  that  equations  (2.5.39)  and 

(2.5.28)  are  identical  if  (y)  ■  {E }  and  if  {£}  -  (0).  Thus  the  {E }  vector  is 
the  stress-displacement  vector  for  the  unforced  vibration  problem,  or  in 
other  words,  { E  (n, <*>,*)}  is  the  depth  dependence  of  the  ntk  fiat  earth  normal 
mode  at  frequency  u>. 

Kazi  (1976)  shows  how  (2.5.39)  can  be  written  as  an  eigenvalue 
problem  for  SH  waves  and  defines  a  Love  wave  operator  whose  eigenvalues 
are  ^K2(d,u)  and  eigenfunctions  are  {^E  (n,u;,x)}  The  definition  of  a  Ray¬ 
leigh  wave  operator  is  not  so  straightforward  because  of  the  P-SV  coupling, 
however,  we  can  still  compute  Rayleigh  wave  eigenvalues  and  eigenfunc¬ 
tions  by  searching  out  the  singular  values  of  (RD(«;,k)].  In  order  to  do  this 
we  will  first  propagate  the  eigenfunctions  from  the  source  depth  to  the  sur¬ 
face  so  that, 

|D(ff  (n,w),0)]  {E  (n,u>,0)}  =  {0}  ,  (2.5.40) 

where 

|TE(A-(n^))] 

[D(JT(n^),0)]  -  -  (2.5.41) 

j®D(  IT  (n,u>),0)j 

Note  that  since  the  { E )  vector  is  a  particular  type  of  stress-displacement 
vector,  it  has  most  of  the  properties  of  the  stress-displacement  vector  and, 
in  particular,  it  can  be  computed  at  different  depths  by  applying  the  correct 
propagator  matrix.  One  property  of  {E }  which  is  not  true  of  (y),  however, 
is  that  {E }  is  continuous  everywhere  with  depth  and  does  not  suffer  a 
discontinuity  at  the  source  depth.  In  fact,  from  equation  (2.5.39)  we  can 
see  that  the  eigenvalues  and  eigenfunctions  are  completely  independent  of 
any  source  characteristics  since  st  in  that  equation  can  be  replaced  by  any 
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other  depth  we  chose  (as  in  equation  (2.5.40))  which  of  course  is  what  we 
would  expect  for  normal  mode  solutions.  It  is  this  decoupling  of  structural 
wave  propagation  characteristics,  as  manifested  by  its  spectra  or  normal 
modes,  from  the  source  and  receiver  characteristics  that  makes  the  normal 
mode  method  an  efficient  solution  of  the  elastic  wave  equation. 

Returning  to  (2.5.40)  we  now  need  to  compute  the  determinant  of 
[D]  which  we  will  refer  to  as  the  characteristic  function 

RA(u,,k)  -  det([RD(u;,k,0)])  ,  (2.5.42) 


-  det([LD(w,k,0)]). 

Assuming  a  free  surface  boundary  condition  we  can  write  the  [D]  matrices 
as  follows: 


and 


(RD(*,M)] 


0  0  10 

M’  M]  M  kH01 


(2.5.43) 


lLD(*,k,0)j 


0  1 
lbd1,0>  lbd,'°I 


(2.5.44) 


We  can  easily  solve  for  the  characteristic  functions  which  are  as  follows: 

RA(w,k)  =  RDjj(u;,k,0)  |f D22(u>,k,0)  (2.5.45) 

LA(w,k)  =  -  ^Dj(u»,k,0)  (2.5.46) 


The  characteristic  functions  are  scalar  functions  of  u  and  k  and  implicitly 
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define  the  eigenvalues  by, 

r^(w*r^  (n>w))  “  ®  *  (2.5.47) 

and 

LA(w,Lif  (n,b/))  «  0  . 

The  normalization  of  the  eigenfunctions  is  arbitrary  and  so  we 
will  assume  a  vertical -displacement  and  SH  shear  displacement  of  unity  at 
the  free  surface.  We  can  solve  for  the  remaining  non-zero  eigenfunction,  the 
P-SV  radial  shear  displacement  in  terms  of  the  [D]  matrix  elements  so  that 

j(n,w,0)  =  1  (2.5.48) 

R£2(n, w,0)  -  -  (jfDa(to/,RJr  (n,to/),0)  R\(w,Rlf  (n,to/),0) 

-  ®D14(to/,Rlif  (n,to/),0)  R 1 D2i(to/,R/f  (n,to/),0)  J 

/(®  D12^>Rjr  M*)  fLDu(“’RK  MM 

■  RDH^’R^ fa’")’0)  R  1 °22^'KK  (wh0)  ) 
j(n,u/,0)  =  0 


R£4(n^i/,0)  “  0 

lE  j(n,to/,0)  =  1  (2.5.49) 

Lr2(n,to/,0)  *  0 

With  the  surface  values  of  the  eigenfunctions  defined,  we  can  compute  the 
eigenfunctions  at  any  other  depth  simply  by  using  the  appropriate  propaga¬ 
tor  matrices. 


In  order  to  use  the  norma)  modes  we  must  say  something  about 
where  the  modes  will  be  located.  As  with  the  branch  points,  for  every 
eigenvalue  at  +k  there  will  be  one  at  -k  due  to  the  fact  that  the  [D]  matrix 
elements  are  functions  of  k*  (or  equivalently  c*).  In  general,  for  all  of  the 
Riemann  sheets  there  will  be  both  pure  real  and  complex  eigenwavenumbers 
(the  exception  to  this  are  SH  and  acoustic  plate  problems  for  which  the 
eigenwavenumbers  are  always  either  purely  real  or  purely  imaginary).  The 
complex  poles  will  be  easy  to  deal  with  since  we  will  include  the  residues  of 
those  complex  poles  which  are  within  the  contour  of  integration,  however 
the  poles  on  the  real  wavenumber  axis  cause  a  problem  since  the  integration 
contour  goes  directly  through  those  poles,  and  we  cannot  know  ofT  hand 
whether  or  not  to  include  their  residue  contributions. 

We  could  compute  the  principle  values  for  these  poles,  but  there 
is  a  simpler  way  to  deal  with  this  problem.  Basically,  we  will  apply  a  per¬ 
turbation  to  the  frequency,  u>,  such  that  the  poles  move  off  of  the  real 
wavenumber  axis  and  can  be  easily  identified  as  being  within  or  outside  of 
the  contour  of  integration.  We  can  allow  the  frequency  to  have  a  small, 
constant  imaginary  component  as  long  as  Im(u/)  <  0  in  order  to  insure  that 
the  Fourier  transform  remains  analytic.  For  each  pole  on  the  real 
wavenumber  axis  for  real  frequency,  we  can  compute  the  group  velocity,  U, 
as  the  slope  of  the  dispersion  curve  in  the  u,  k  space,  or, 


U(n,ur)  *  da//dk 


w,RJf  (n,u >)’ 


(2.5.50) 


A  small  change  in  w,  6u,  will  thus  cause  a  small  change  in  rA(d,v),  SK , 
such  that, 


SK  =  iu/ U  . 


(2.5.51) 


So  for  poles  on  the  positive  real  wavenumber  axis,  if  the  group  velocity  is 
positive  then  a  small  negative  imaginary  perturbation  of  frequency  will 
move  the  poles  into  the  fourth  quadrant  where  they  would  be  within  the 
integration  contour,  and  if  the  group  velocity  is  negative  then  the  frequency 
perturbation  will  move  the  poles  into  the  first  quadrant  where  they  will  be 
outside  of  the  integration  contour.  However,  for  every  pole  at  +k  there  will 
be  one  at  -k  and  it  is  easy  to  6how  that  a  pole  at  +k  with  group  velocity  U 
will  have  a  companion  pole  at  ~k  with  group  velocity  -U.  So  the  companion 
poles  to  those  on  the  positive  real  wavenumber  axis  with  positive  group 
velocities  will  have  negative  group  velocities  and  will  move  into  the  second 
quadrant  where  they  will  be  outside  of  the  contour  of  integration.  The 
companion  poles  to  those  on  the  positive  real  wavenumber  axis  with  nega¬ 
tive  group  velocities  will  have  positive  group  velocities  and  will  move  into 
the  third  quadrant  where  they  will  be  within  the  contour  of  integration. 
The  net  result  is  that  all  poles  with  positive  real  wavenumbers  and  positive 
group  velocity  will  contribute  their  residues  to  the  wavenumber  integral  and 
those  poles  with  positive  real  wavenumbers  and  negative  group  velocity  will 
contribute  with  their  companion  poles  at  -k.  The  wavenumber  integration 
contour  in  the  complex  wavenumber  plane  along  with  the  branch  points, 
branch  cuts  and  poles  are  shown  in  figure  2-2  for  the  genera]  P-SV  half 
space  problem.  We  can  see  from  figure  2-2  that 

+00 

/  +  Ra1  +  Kfi1  =  -  2*  i  residues,  (2.5.52) 

JO  D 

where  the  arcs  at  infinity  do  not  contribute  and  the  sum  of  residues  are 
those  within  the  integration  contour,  f. 
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Figure  2-2.  Wavenumber  integration  contour  for  a  complex  frequency  with 
a  email,  negative  imaginary  component. 


We  now  turn  our  attention  to  the  evaluation  of  the  residues. 
Returning  to  equation  (2.5.28),  we  can  write  the  solution  for  the  stress- 
displacement  vector  at  the  top  of  the  structure  as  we  did  with  the  eigen¬ 
functions 


(D(0)]  {y(0)> 


{0} 

|BD(0)|  |a(0j,):  <S) 

i 


(2.5.53) 


Once  gain  we  assume  a  free  surface  so  that  (D(0)]  is  given  by  equations 
(2.5.43)  and  (2.5.44).  Substituting  these  relations  in  (2.5.53)  it  follows  that, 


Ryi(°) 

l 

Ry2(o) 

R* 

RB»1,<0) 

-rDs1(o) 


-rd.2(°> 


lRBD(0)i  IrMO.i,)]  {„£}  12.5.54) 


Ry5(°)  *  Ry4(°)  *  0 » 

and 

Ly,(°) - ^  l®D(0);  |tA(0^)]  (LE) ,  (2.5.55) 


Ly2(°)  -  0. 

We  can  see  that  the  two  by  two  minors  of  the  two  by  four  P-SV 
(R®D]  matrix  appear  repeatedly  throughout  the  analytical  development  of 
the  eigenvalues,  eigenfunctions  and  the  stress-displacement  vector.  In  order 
to  save  writing  we  will  define  the  four  by  four  anti-symmetric  minor  matrix. 
[M],  as, 

Mij(w,k,i)  -  RDjj(w,k^)  RD2j(u;,k,i)  (2.5.56) 

-  RDlj(w’M  RD2i(w,k,r)  • 


Obviously, 


69 


M ifivM  •  -  M-fu/^,*) , 

MII  "  J  “  M«  "  M44  “  °- 

We  will  also  have  two  versions  of  the  minor  matrix,  [BM]  and  [Tm] 
corresponding  to  [^D]  and  [j^D].  The  minor  matrix  can  now  be  used  to 
simplify  the  solution  to  equation  (2.5.54)  as  follows, 

R»,(0)  I  pMutO)  0  -  "m^O)  -  BMjJO)] 

*y»(°>[  '  R*  |  0  Bm1i(0)  Bmij(0)  bMu(0)  j  IrA(0^1  <re>  ■ 

(2.5.57) 

We  can  also  see  from  equation  (2.5.45)  that, 

„A  -  BMjjfO) .  (2.5.56) 

The  P-SV  radial  displacement  eigenfunction  at  the  surface,  given  by  equa¬ 
tion  (2.5.48),  can  be  written  as  follows 

K’M  -  2-r; - m  -  (2.5.59) 


BMu(u>,RiSf(D^),0) 
Bmu(u*RX‘  (n^),0) 


bM18(w,rJT(d^),0) 
BM„(w,Rff  (n,»),0)  ’ 


where  Re(n,w)  is  the  Rayleigh  wave  surface  ellipticity  for  the  nll>  mode  at 
frequency  u. 

Poles  of  equation  (2.5.52)  will  occur  whenever  RA  ~  0  which 
implicitly  defines  the  eigenwavenumbers.  We  can  simplify  equation  (2.5.57) 


eigenfunctions,  but  first  we  must  state  two  genera)  properties. 


Mu(w,k4)  -  -  M24(^Jc4)  (2.5.60) 


This  is  true  for  all  w,  k  and  i.  It  is  obviously  true  for  the  boundary  condi¬ 
tions  given  by  equations  (2.5.16)  through  (2.5.20)  at  the  top  and  bottom  of 
the  structure,  and,  with  much  tedious  algebraic  manipulation,  it  can  be 
shown  to  be  true  at  all  depths.  The  second  property  allows  one  to  relate 
the  elements  of  a  P-SV  propagator  matrix  which  propagates  the  stress- 
displacement  vector  upward  between  two  depths  to  the  elements  of  the 
downward  propagator  matrix. 


!rA(*,*0)J 


RASs(*0’*)  RA4s(,0’*i  “  RAll(*0*»>  “  RA2s(,0’1) 
rAj4(*o»*)  RA4<(*0’*)  “  RA14(*0’*)  “  RA24(*0’*) 

“  RA3l(*0’*)  ”  RA4l(*0»*)  RAll(*0’*^  RA2l(I0’*) 

~  RA22(*01*^  ~  RA42^,0’*^  RA12(*0’*^  RA22(,0**) 

(2.5.61) 


This  can  be  shown  by  using  the  relations  that  exist  among  the  elements  of 
the  P-SV  layer  propagator  matrix  (equations  (2.4.44))  and  then  reversing 


the  order  of  layer  multiplications.  Equation  (2.5.6)  is  a  general  property  of 
the  propagator  matrix;  it  is  true  for  all  u>,  k,  z,  and  zQ  and  it  is  true  even  for 
an  arbitrarily  inhomogeneous  structure  with  depth.  We  now  define  the 


numerator  vector  in  (2.5.57)  as 


Ryi(*o) 

i  , 

' 

rni« 

Ry2(*o) 

=  ra 

RN2(«) 

, 

(2.5.62) 


Using  equations  (2.5.58),  (2.5.59),  and  (2.5.60)  we  can  solve  for  [N] 
evaluated  at  an  eigenwavenumber  as, 
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r«.<0' 

rN,<#> 


l-  **(•*•) 


00  1,. 

0  0  R'  R‘* 


IpA(0^,)l  {*£) 


k-Rlf  (n^) 
(2.5.63) 


Using  (2.5.61)  this  can  be  expressed  as  follows, 


rN  ,«°» 


-  ®M,,(0)  11  H«1 


(2.5.64) 


RASl(*,°)  -  RA41(«,0)  RAn(«g  0)  RA21(*,0) 
“  RA42(**,°)  RA12^V0^  RA22(*«,°) 


•<R*> 


rN,'01 


k-Rlf  (■,«) 


k-Rlf(n^) 


Re 


k-RA(n,w) 


We  can  write  the  Rayleigh  wave  eigenfunctions  at  the  source  depth  as 


(r£  (n,u>,*t)} 


ra1|(V°)  RA12^*i’°) 
RA2l(*»’0)  RA22(V°) 
RASl(*»,°)  RA32^*»’°) 

RA4l(*,»°)  RA42^*»’°) 


i  1 

k-R*(n,u/)  iR*^) 


(2.5.65) 


If  we  redefine  a  new  vector  for  the  source  jump  vector  as, 


Ir£(o,w)]  *  [r£j  rE4  -rEj  -rE2] 


k-R/f(n,ui) 


(2.5.66) 


and  using  (2.5.65)  we  can  write  (2.5.64)  as 
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rN,|0) 

rNS'°> 


.n<°I 

Rn4 


(2.5.67) 


k-R/f(«,w) 

=  -  BM„(0)  k.  ([RS(n,w)]  (RE(n,u;,iJ}){RE(n,a;,0)} 


where  we  have  expanded  {RN}  to  four  P-SV  components.  Finally,  we  can 
write  the  numerator  vector  at  any  other  depth  by  applying  the  propagator 
matrix  so  that. 


(rN(«)> 


k-RJf  (n,w) 


(2.5.68) 


-  -  BM„(0) 


([R£(n,u>)j  (RE(n,uMf)}  ){RE(n,w,*)} 
k -R#f  (n,w)  '  ‘ 


We  can  also  define  a  numerator  vector  for  Love  waves  and  without  repeat¬ 
ing  the  detailed  derivation  one  can  show  that, 


{iN(*)> 


k «!*  (n,w) 


(2.5.69) 


-  -  BDj(0) 


k«Lff  (n,w) 


([L2(n,w)]  (LE(n^^#)>  J{LE(n,w,i)} 


where 


{LE(n,w)]  =  [LS2  -  lE,] 


k-LK(D,w) 


(2.5.70) 


and 


{i,y(«)>  -  -a-  • 


(2.5.71) 


r 


i 
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We  may  now  evaluate  the  wavenumber  integrals  given  by  equa¬ 
tions  (2.5.36)  and  (2.5.37)  in  terms  of  branch  cut  integrals  and  residue  con¬ 
tributions.  Using  equations  (2.5.38),  (2.5.52),  (2.5.58),  (2.5.68),  (2.5.60)  and 
(2.5.71)  we  can  write  the  frequency  dependent  displacements  at  a  receiver 
location  (rf,  $T,  *r)  due  to  a  source  confined  to  a  horizontal  plane  at  depth  *# 
as  follows, 

Ru(w^r,8r4r)  =  -  RaI  “  R^I  (2.5.72) 


“‘HE  (rA  (n,u;)  |R2(n,w,m)]  {RE(n,w,sf)}R*(n,w,m,rr,8r,*r)J  , 


and 


Lu(u>,rr,0r,ir)  =  -  l0l 


(2.5.73) 


“  >E  E  (lA(“’w)  kE  L*  (n^»m>rr^r’Ir))  * 


where  the  A ’s  are  scalar  amplitude  factors  and  are, 

k  “MjjIO) 


,  A  (n,w)  -  - 


0RA/dk 


(n,w) 


(2.5.74) 


and 


LA(n,w)  =  - 


k  “•■DjIO) 

aLA/ak 


k-LK  (n,w) 


(2.5.75) 


rA  and  lA  are  defmed  by  equations  (2.5.45)  and  (2.5.46),  [RE]  and  [LE]  are 
defined  by  equations  (2.5.66)  and  (2.5.70),  and  R$  and  are  defined  as 
follows, 


i 


(2.5.76) 


R*(n,u>,fn,rr#r,tr)  -  Rr,(n,u/,sr)  ^(n,u;(in^r,#r) 

+  R£2(n,w,ir)  6(n,w,m^r,^r)  , 

l4>  (n,w,m,rr  tfp4r)  -  ,(n,u;,ir)  fc(n,w,m,rr,0r)  ,  (2.5.77) 

where  P,  6  and  C  are  modified  vector  cylindrical  harmonics  and  are, 


P (n,w,m,rr,tfr)  =  e,  (krp)  e1”1^ 


k-R#f  (n,w) 


(2.5.78) 


6(n,ta»,m,rr,tfr)  *  ef 


aHm2^  (kr)  im*r 


d(lcr) 


k«RJf  (n,to»),  r-rf 


(2.5.79) 


+ 


Hm2)(krr)  btm8 


(kr)  be 


k«RJf  (n,w),  9~9 T 


and 


C(n,u>,m,rr,0r)  -  ep 


<>  (k',1  a,™* 


(krr) 


be 


k-RK  (n,w),  9~9f 


(2.5.80) 


aw^l  (kr)  in,,, 

d(k0  * 


|k-=LA'(n,w),  r-rr 


2.6  The  Branch  Cut  Integral  Contributions 

Finally  we  turn  our  attention  to  the  branch  cut  integrals  RaI, 
R^I,  and  L*1-  Let  us  first  consider  the  Rayleigh  wave  branch  cut  integrals 
as  shown  in  figure  2-2.  In  order  to  facilitate  the  evaluation  of  the  branch 
cut  integrals  we  define  the  following  complex  valued  functions  of  the  real 
positive  scalar  variable  17. 
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**  P-wave  branch  cut,  (2.6.1) 

(ti,u)  *  S- wave  branch  cut, 

where 


RqK  (0,oj)  -  P-wave  branch  point, 

-  W/a<N> 

R0*  (0,u>)  =  S— wave  branch  point, 

=  u,/*<N> 

and  0  <  t)  $  oo  . 

We  will  also  denote  wavenumber  values  immediately  to  the  right  and  left  of 
the  branch  cuts  as  viewed  in  figure  2-2  with  +  and  -  superscripts.  We  can 
now  write  the  branch  cut  integrals  as  follows, 

o 

R.1  *  /  «R„*  kJ  '  <*1  (2  6  2) 

00 

00 

+  f  f(Ratk~('7’ta'))  Ro* ' d*  ’ 
and  similarly  for  the  g^I,  where 

and  f(k)  is  the  wavenumber  integrand  function.  We  can  combine  the  two 

integrals  in  (2.6.2)  to  obtain  the  following 

•  • 

00 

-  R.1  -  |  («R.*  +(->.“))  -  f<R„*  '(».")))  r„K  -il  (2.6.3) 
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The  P  and  S-wave  branch  cuts  define  discontinuities  in  the  bot¬ 
tom  half  space  vertical  wavenumbers,  and  so  that  as  one  crosses 
the  P-wave  branch  cut,  i/jN)  -*  -  and  as  one  crosses  the  S-wave 

branch  cut,  i/jN)  -*  -  The  only  factors  in  the  wavenumber  integrand 

functions  which  will  be  discontinuous  across  the  branch  cuts  will  be  the  ele¬ 
ments  of  the  stress-displacement  vector,  and  so  in  order  to  evaluate  (2.5.29) 
we  need  to  compute  the  stress-displacement  difference  function  across  the 
branch  cut  which  we  define  as  follows, 

<^y(*7,*)>  -  {y+(*»,*))  -  {y“(*7**)> ,  (2.6.4) 

where 

{y+(«M)}  *  {Ry(M)> 

k-R*+(9) 

and 

{y"  (»?,*)}  =  (Ry(M)} 

k-R  K-(„) 

For  Rayleigh  waves,  we  can  use  equation  (2.5.57)  to  compute  (Ry  +  )  and 
{Rv-}  at  the  surface.  Remembering  that  the  propagator  matrix  and  source 
jump  vector  are  both  single  valued  functions,  we  can  compute  the  stress- 
displacc  :nent  branch  cut  jump  at  the  surface  as, 

.  L J-*se»l  J.!^l 

<Ry,(l.«)  ”  ( ®M,,(0)  |  I  ®M..(0)  I 

)  o  0  <  i  14 - 

bm12(°)  ]  ®M„(0) 

- 


^(O)  J  *1^(0)  ^(O) 

From  equation  (2.5.60)  it  is  easy  to  show  that 


B 

M„(0) 

B 

i 

M„(0) 

We  can  see  that  equation  (2.6.5)  resembles  equation  (2.5.63) 
which  expresses  the  solution  of  the  stress-displacement  vector  in  terms  of 
the  discrete  spectra  of  the  Rayleigh  wave  operator.  In  order  to  derive  a 
similar  solution  for  the  branch  cut  integrals,  we  need  to  determine  the 
improper  eigenfunctions  which  will  constitute  the  continuous  spectra  of  the 
Rayleigh  wave  operator.  We  will  do  this  by  first  considering  the  stress- 
displacement  vector  (Rx)  which  obeys  the  radiation  condition  in  the  bottom 
half  space,  but  not  necessarily  the  free  surface  boundary  condition.  In  com¬ 
puting  the  (Rx)  vector  we  will  also  assume  no  sources  so  that  from  equation 
(2.5.26)  we  can  write  the  following. 


rdii  rdi»1  |r*i 


RUM  Ru27 


R  D1S  R  D14  R*S 
R  DM  R  D24  R*4 


(2.6.7) 


Equation  (2.6.7)  will  hold  at  all  depths  and  so  we  will  evaluate  Rx2  and  Rx2 
at  the  surface  where,  using  the  (D)  matrix  minors  we  arrive  at, 


R*,(°) 

R*2<°) 


-  BMjj(0)  -  ^ 


®M„(0)  DM1S(0) 


V°) 

14<°) 


RXj(°) 

R*4(°) 


(2.6.8) 
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Of  course,  equation  (2.8.8)  is  only  true  if  bMm(0)  *  0.  We  will  now  assume 
that  Rxs(0)  -  0  so  that  we  can  solve  for  rx^O)  and  rXj{°)  of 

RX4(0).  We  will  also  evaluate  these  at  Rff  +(*)  and  RX  ~{v)  so  that. 


**>+(0)  •  ^ »*• (0) ' 

R  *  bm^(0)  k 

and  similarly  for  R*f(°)  and  RXj"(0).  We  will  specify  that, 

rx4+(o)  *  rx4-(°)  ’ 


(2.6.9) 


(2.6.10) 


so  that, 


and 


BM2‘t(0)  BM~(0) 

rXi  (0)  =  bm+(o)  bMj"(o)  rX* 


®M.+(0)  BM,-2(0)  _ 

rx2+(0)  -  ~rr  •  r**  <°> 


(2.6.11) 


®M+(o)  bm-(0) 

We  can  now  define  the  stress-displacement  vector  {R^}  as  follows, 

PrX)  _  fox*)  ~  <RX~) 

(rE)*  <rx,(0)  ■  {Rv,+(0»  -  iRxr«»>  ’ 


(2.6.12) 


so  that 


R*l<°> 

1 

RX2+(°)  ~  R*2  W 

R*7<°) 

< 

rX!+(0)  -  RXf(O) 

R*#> 

R*4«>) 

0 

0 
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G  . 


where 


"M#) 

bmm(0) 

bm34(o) 

bm..(0) 


(2.6.13) 


Following  the  same  procedure  except  setting  Rx4+(0)  -  R*4  (0)  ■  0  we  can 


show  that 


BM38(0) 

®M„(0) 


(2.6.14) 


Equations  (2.6.13)  and  (2.6.14)  along  with  (2.6.6)  allow  us  to 
express  equation  (2.5.5)  in  the  following  manner. 


(2.6.15) 


Ryi(°)|  (  -  BMm(0)  |  f°  0  1  Rf* 


Ry*(0)|  ^,,(0)  0  0  Re  rc! 


IrA(<M.)] 


We  can  now  repeat  the  derivation  of  equations  (2.5.63)  through  (2.5.68)  to 


show  the  following. 


(2.6.16) 


(*Ry(»M)}  *  s  -  |(r£(v^)]  {r*<,.,.>)|  |r^  (*»".*)  J  » 
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where 


and 


Ir2(^)1  "  Ir£j  rs4  -r£i  ”r£s) 


k-R  k[n*>) 


(2.6.17) 


{r£(»?  .".*)}  *  IrA(*»°)] 


k“R^(n.w) 


(r£  (»?,w,0)}. 


(2.6.18) 


In  order  to  clean  up  the  notation,  we  have  dropped  the  a  and  ft  subscripts 
denoting  P  and  S-wave  branch  cuts  throughout  the  developments  but  it  is 
understood  that  there  will  be  two  versions  of  (2.6.15)  corresponding  to  the 
two  branch  cut  integrals,  RaI  and  R^I.  Following  a  similar  analysis  for  the 
Love  waves  we  can  express  the  stress-displacement  vector  jump  across  the 
Love  wave  branch  cut  as  follows. 


tfLy(»?i*)} 

where 


(2.6.19) 


[L£(f*,w)J  =  |le2  -lEjJ 


k-Ltf  (*,*)’ 


(2.6.20) 


and 


(i?,a),*))  -  1lA(*,0)] 


{lE(*,u;,0)>  , 

k-I *(■»,*) 


(2.6.21) 


{l£(»;,w,0)) 
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